The importance of DNA sequence for nucleosome positioning in transcriptional regulation

Abstract

Nucleosome positioning is a key factor for transcriptional regulation. Nucleosomes regulate the dynamic accessibility of chromatin and interact with the transcription machinery at every stage. Influences to steer nucleosome positioning are diverse, and the according importance of the DNA sequence in contrast to active chromatin remodeling has been the subject of long discussion. In this study, we evaluate the functional role of DNA sequence for all major elements along the process of transcription. We developed a random forest classifier based on local DNA structure that assesses the sequence-intrinsic support for nucleosome positioning. On this basis, we created a simple data resource that we applied genome-wide to the human genome. In our comprehensive analysis, we found a special role of DNA in mediating the competition of nucleosomes with cis-regulatory elements, in enabling steady transcription, for positioning of stable nucleosomes in exons, and for repelling nucleosomes during transcription termination. In contrast, we relate these findings to concurrent processes that generate strongly positioned nucleosomes in vivo that are not mediated by sequence, such as energy-dependent remodeling of chromatin.

Introduction

The nucleosome is a foundational protein complex in eukaryotic organisms formed by an octamer of pairs of histone subunits (two pairs of H2A/H2B and a tetramer of H3/H4). This complex wraps 147 base pairs (bp) of DNA (1, 2). The positioning of nucleosomes is very important for the regulation of transcription, and the relevance of nucleosome positioning for the regulation of genes can be viewed from many perspectives. The broadest unit to span interaction potential are topologically associated domains (TADs). These are characterized by a loop compartment that brings together potentially interacting proteins at enhancers and promoters of genes close to each other (3). The boundaries of these logical units are set up at so-called insulator regions, which are defined by an interplay of CTCF binding (4) and a consequential strict positioning of a well-positioned nucleosome array (5). Whether the genes inside such units are transcribed is an intricately regulated process involving nucleosome positioning at all levels. The binding of transcription factors (TFs) to promoter and enhancer regions is crucial for the initiation of transcription. As a means of controlling access to transcription factor binding sites (TFBSs), nucleosomes occupy DNA and need active displacement through (pioneer) TFs or energy-dependent remodeling to avoid random triggering of gene expression (6, 7). To achieve an orderly rate of transcription, the RNA polymerase II (RNAPII) is highly regulated by nucleosomes. In particular, the +1 nucleosome directly downstream of the transcription start site (TSS) and a sequential array of well-positioned nucleosomes are needed to reliably initiate transcription (8, 9, 10). Along its way, RNAPII encounters more fuzzily positioned nucleosomes in intronic regions and well-positioned nucleosomes in exons.

Because exonic DNA is the template for the subsequent translation of proteins, the positioning of nucleosomes supports the faultless transcription of these regions by serving as roadblocks to decrease the speed of RNAPII elongation, by marking alternative splice sites (11), and by triggering feedback signals for other concurrent transcription initiation events on the same gene (12).

The transcription ends at the transcription termination site (TTS), which sets up the signal where polyadenylation will be applied to the resulting pre-mRNA. These polyadenylation sites (poly[A] sites) are known to exhibit a considerate nucleosome depletion directly around the DNA-encoded prospective polyadenylation site compared with their nucleosome-occupied surroundings (13).

An extra layer of nucleosome functionality in all of these stages is added through the epigenetic modifications of histone proteins. Depending on the specific biochemical modulations of the histone subunits, nucleosomes can carry marks for the active recruitment of further proteins and show a different flexibility in terms of the energy needed for displacement (14, 15). Because the differential increase or decrease in chromatin accessibility at all of these specific functional sites regulates transcription and thus cell identity, it is crucial to understand the principles that determine nucleosome positioning.

The influences on the positioning of nucleosomes have been debated for many years and can be grouped roughly into DNA sequence, trans-acting factors, and active chromatin remodeling (16, 17, 18). The role of the DNA sequence is particularly controversial. It was postulated that there is a periodic repeat of A/T dinucleotides that support the binding of nucleosomes (19, 20, 21). This pattern has been extended to be represented by local DNA structure (22) with the help of DNAshape (23). The benefits of the supposed patterns for nucleosome support are linked to the flexibility of the DNA strand to facilitate an easier attraction of nucleosomes (24, 25), which can be seen particularly well in extremely stiff poly(dA:dT) tracts that counteract nucleosome binding (26).

Some of these findings have been discovered with the help of a variety of computational tools to predict nucleosome positions, each with a different focus of training data, prediction target, and methodology (27, 28, 29, 30). Among these, there are tools that employ HMMs (31, 32), SVMs (33), and recently more and more convolutional neural network (CNN) architectures (34, 35, 36, 37). The difference in meaning between sequence-intrinsic signals and the actual in vivo positioning of nucleosomes is fundamental to understand. The sequence can only explain the potential of any given position to attract nucleosomes. Whether it is actually occupied by nucleosomes or opened up by regulatory mechanisms cannot be determined without the help of experimental chromatin accessibility data for the particular cell type. In other words, although experimental data can provide insights into the actual nucleosome positioning within cells, it does not directly elucidate the underlying reasons or mechanisms behind such positioning. On the contrary, machine learning on DNA sequence can give an intrinsic estimate of nucleosome attraction potential but cannot answer the question where nucleosomes are located in vivo. Therefore, we need to integrate both information to estimate the overlap of sequence-intrinsic nucleosome support and experimental chromatin accessibility evaluation (38). For that purpose, we use nucleosome occupancy data (MNase-seq) (39), chromatin accessibility data (DNase-seq) (40), and ChIP-seq data for both histone modifications and TFs, mainly in the two cell lines GM12878 and K562 available from ENCODE (41).

In this study, we evaluate the importance of DNA sequence for nucleosome positioning at all elements involved in transcription in the human genome. For that purpose, we developed the nucleosome formation score (NF score), which describes the sequence-intrinsic support for attracting nucleosomes purely based on DNA. The score is derived using a classification algorithm, which is applied in a genome-wide manner. We provide it as a simple data resource, which can be used without further computation for similar analyses in the future. Because of its relative robustness against overfitting, we created a random forest (42) classifier based on the Guo2014 dataset (33) of nucleosomal sequences, which is derived from MNase-seq in human CD4⁺ T cells (39). Our classifier uses a transformation from raw DNA sequences to a frequency spectrum of local DNAshape (propeller twist, helix twist, minor groove width, electrostatic potential, obtained with the DNAshapeR package (43)). This binary classifier was applied in a sliding-window manner to the human genome in two different resolutions, yielding a genome-wide mapping of sequence-intrinsic nucleosome support.

We aspire to give a comprehensive overview of the crucial parts of the transcriptional machinery and provide our genome-wide measure for sequence-intrinsic nucleosome support as a simple data resource.

Our findings are embedded into a framework of previous studies that showed the relevance of DNA for nucleosome positioning for some specific elements, and we add some new insights about the sequence definition of nucleosome positioning in transcription. We demonstrate that there is a clear influence of DNA sequence for supporting nucleosome binding at positions of high competition with other DNA-binding proteins and where it is beneficial to accumulate nucleosomes by default, such as exons, or repel them in TTSs. In contrast, we give examples of positioned, but sequence-independent, nucleosome arrays in proximity to the TSS or around insulator sites and put our findings into the context of more unspecific, fuzzy nucleosomes.

Results

Sequence-intrinsic nucleosome support characterizes functional genomic regions

We constructed a resource to estimate the support of DNA sequence for attracting nucleosomes based on a random forest classifier trained on nucleosome occupancy data. This nucleosome formation (NF) score was applied to a range of different genome locations. A broad range of such entities can be found in the histone ChIP-seq–derived mappings of chromHMM for the cell lines K562 and GM12878 (44). Fig 1A shows the mean sequence-intrinsic nucleosome support (NF), the nucleosome occupancy (MNase-seq), and the chromatin accessibility (DNase-seq) for the 147 bp in the center of any given chromHMM fragment. In total, there are 622,257 fragments in K562 and 571,339 in GM12878. A list of the individual fragment numbers for each subgroup can be found in Table S1. We grouped all chromHMM states into the three categories: “regulatory” (promoters/enhancers), “heterochromatin/repetitive,” and “transcription” (transcriptional transition/elongation, etc. and insulators).

Nucleosome support in different regions of the genome and around histone modification ChIP-seq peaks.

(A) chromHMM. Nucleosome support (NF score), nucleosome occupancy (MNase-seq), and chromatin accessibility (DNase-seq) in all states of the chromHMM mappings. The highlighted regions exhibit a high NF score in active regulatory regions (active and poised promoters) and in insulators, as opposed to low nucleosome support regions in heterochromatin. Active promoters are rather accessible, and poised promoters are closed off mostly. Heterochromatin exhibits pronounced inaccessibility, whereas insulators demonstrate hybrid accessibility characterized by relatively elevated levels of both DNase-seq and MNase-seq signals. (B) Nucleosome support of histone modifications. A distinct peak is observed in proximity to each histone modification ChIP-seq peak, providing evidence of the influence of DNA sequence on the positioning of individual nucleosomes. This influence is higher in marks indicating a regulatory function such as poised promoters (H3K27me3, green), active promoters (H3K4me3, blue), and more generally accessible chromatin (H3K27ac, red) as compared to repressing histone marks (H3K9me3, purple).

Table S1

The most intriguing feature is the relation between sequence-intrinsic nucleosome support and actual nucleosome occupancy. Active promoters show the highest DNA-intrinsic nucleosome support (NF: ∼0.8). However, being crucial places for transcriptional regulation, they are also the most accessible regions (low MNase-seq of 0.4–0.55 and high DNase-seq of 0.65–1). The same nucleosome support can be observed in poised promoters, but with closed chromatin (low DNase-seq, high MNase-seq). This observation can be put in contrast to the highlighted group of heterochromatin fragments. Here, we observe rather low nucleosome support in both cell lines of roughly 0.4. At the same time, the locations are very inaccessible with a nucleosome occupancy of 1 and a very low DNase-seq signal of 0.03.

The last highlighted group in Fig 1A are insulators. In this group, we can observe a relatively high nucleosome support of 0.6–0.66 and a medium accessibility with relatively high values for both MNase-seq of 1 and DNase-seq signal of 0.2–0.3.

As a side note, it can be observed that enhancers show the most variability in their status. This effect on the NF score is most prominent in the groups of Weak_Enhancer (groups 6 and 7), in which the NF score varies almost about 0.2 between the minimum in GM12878 and the maximum in K562. These groups are the most ambiguously defined ones, however, based on their histone marks (44). Nevertheless, also in the Strong_Enhancer group (group 4) and the Active_Promoter group (group 1), there is a comparatively high difference in chromatin accessibility. This could be due to experimental biases but because the direction of the differences is not the same in both instances.

In summary, these observations suggest that regions involved in active regulation of transcription have significantly higher levels of nucleosome support with well-defined sequence context. Whether these important locations are actually occupied by nucleosomes is dependent on whether they are used in vivo (active versus poised promoters). By comparison, constitutive heterochromatin has a much lower sequence-intrinsic nucleosome support while being mainly inaccessible.

Histones can carry different kinds of modifications that determine their function and the flexibility of the surrounding chromatin accessibility (14, 15). Because the chromHMM mapping is derived from ChIP-seq data of histone modifications, we can gain more detailed insights into the local nucleosome support relative to the actual position of the histone modification. Therefore, we analyzed 518,609 histone ChIP-seq peaks from ENCODE for the GM12878 and K562 cell lines. We chose the histone marks according to their role in relation to the derived chromHMM states. The H3K4me3 peaks are representative of active promoter histone marks (GM12878: 194,130; K562: 240,828), H3K27ac is representative of enhancers (GM12878: 114,738; K562: 102,352), H3K27me3 represents poised promoters (GM12878: 100,266; K562: 477,102), and H3K9me3 is repressive chromatin as a general indication for constitutive heterochromatin (GM12878: 37,756; K562: 11,778). Whenever there are multiple equivalent experiments for the same target, we calculated the mean over all experiments in a single line.

Fig 1B shows the average NF score around the set of different histone modification ChIP-seq peaks. It can be seen that, in the specific experimental configurations, the previous findings of a higher NF score in regulatory important regions can be confirmed. For both cell lines, it can be observed that the three activating marks are higher than the heterochromatin mark. In addition, all histone modifications show a clear peak of nucleosome support around the central ChIP-seq point-source peak. The poised promoter mark H3K27me3 is consistently the most sequence-positioned in both cell lines and can be found at NF score levels of 0.82 and 0.86, respectively. The next highest pattern on both sides is the active promoter mark H3K4me3, which shows a peak at an NF score of 0.67–0.78. The general accessibility mark in promoters and enhancers H3K27ac is next at 0.64–0.76, and the peak of the heterochromatin mark H3K9me3 is located at the lowest NF score level of 0.44–0.66. These results indicate a local sequence preference for nucleosome binding around the central peak position of the histone modifications with a distinctly higher nucleosome support for regulatory elements. Because the modified histones are part of the nucleosome, it can indeed be expected that the NF score is highest around its central binding places, but it is remarkable that the sequence-intrinsic influence is highest at putative positions of interaction between nucleosomes and proteins such as TFs.

Regulatory nucleosomes: competition between sequence-intrinsic nucleosomes and transcription factors

Promoter regions

The previous analyses suggest that DNA structure plays different roles in nucleosome binding in different biological environments. To enhance our comprehension of these distinctions, we direct our attention toward regulatory loci, specifically encompassing all human promoters. This analytical approach allows us to establish correlations between the unchanging attributes of sequence-intrinsic regulatory potential and its consequential influence on gene regulation. The promoters are defined as 1 kb upstream of all 30,142 RefSeq-defined human genes (see the Materials and Methods section). Alternative TSS locations are counted as individual promoters, and strand direction is aligned to 5′-3′ direction.

We applied the high-resolution NF score (see the Materials and Methods section) to all human promoters and compared this cell-type unspecific signal with the actual experimental evidence of chromatin accessibility in K562 and GM12878, measured by ENCODE MNase-seq, DNase-seq, and H3K4me3 ChIP-seq data, as well as a derived DHS score, which integrates DNase-seq peaks over 403 primary cell lines (see the Materials and Methods section). In addition, we include PhyloP, a bp resolution evolutionary conservation score.

Fig 2A shows the position-specific mean distribution for all available genes in relation to the TSS. An extended version of this figure showing 1 kb up- and downstream of the TSS can be found in Fig S1. The overall trend of the different data sources is very similar between both cell lines and confirms the general expectations of chromatin accessibility. Although the general occupation with nucleosomes declines toward the TSS (MNase-seq from 0.9 to 0.5), the amount of promoter histone modification H3K4me3 increases toward the TSS (from 10 to 25 in K562, from 4 to 7 in GM12878) with a distinct valley about 100 bp upstream of the TSS. This decline in nucleosome occupancy can be explained simultaneously with the gain of chromatin accessibility by putative proteins such as TFs, which is described by a peaked rise of DNase-seq (from 0.25/0.3 to 1.0/1.6) and a DHS score (from 10 to 145) in the same area. The NF score increases toward the TSS with a more steady, less peaked slope from ∼0.5 at −1 kb to 0.85 at the TSS. The evolutionary conservation (PhyloP) stays relatively stable at a level of 0.09 until it increases steeply to 0.45, starting at −200 bp and continuing with a high evolutionary conservation downstream into the gene.

Chromatin accessibility of human promoters.

(A, B) Profiles for experimental chromatin accessibility (MNase-seq, DNase-seq, H3K4me3 ChIP-seq) and static data (green: DHS score, NF score, PhyloP) as (A) the mean of all human RefSeq promoters and (B) the promoters of the individual genes NANOG and TP53 in K562 (blue) and GM12878 (red). The dotted line at position 0 marks the TSS. (A) Overall promoter regions. The chromatin accessibility increases, whereas the general nucleosome occupancy decreases. The sequence-intrinsic nucleosome support (NF score) increases toward the TSS along the increasing chromatin accessibility. Furthermore, the H3K4me3 ChIP-seq and MNase-seq signals exhibit an NDR upstream of the TSS, and the evolutionary sequence conservation represented as PhyloP score is generally decreased in the promoter region but increases steeply toward the TSS. (B) NANOG profiles. In the context of the two cell lines, NANOG, because of its role as a developmental pioneer factor, does not exert any influence in K562 and GM12878. The marked area denotes an overlap of high sequence-intrinsic nucleosome support and a high MNase-seq peak. That means the main site of nucleosome attraction potential is occupied by well-positioned nucleosomes. TP53 profiles. TP53 is an important tumor suppressor in both cell lines. The NF score follows the peaks of high accessibility in the DHS score and DNase-seq in location. In the annotated instance, the region characterized by elevated nucleosome support exhibits a nucleosome-free state in vivo, indicating a displaced nucleosome position. Both gene promoters show an increased PhyloP score at the high NF score sites, indicating an evolutionary conservation of the nucleosome-supporting sequence. Overall, it should be noted that the NF score correlates positively with the accessible chromatin data (DHS, DNase-seq) in actively used promoters but correlates with nucleosome positioning data (MNase-seq) in closed promoters, demonstrating the original nucleosome positioning potential by the DNA sequence, which is overwritten upon promoter activation. y-Axes: MNase—density graph of signal enrichment; H3K4me3—control normalized tag density; DNase—read depth normalized signal; DHS—DHS score (see the Materials and Methods section); NF—NF score (high resolution); PhyloP—PhyloP score.

Promoter profiles.

Profiles for experimental (MNase-seq, DNase-seq, H3K4me3 ChIP-seq) and global (green: DHS score, NF score, PhyloP) data as a mean of all human RefSeq promoters and the promoters of the genes NANOG and TP53 in K562 (blue) and GM12878 (red). This plot is an extension of Fig 3 in the main article, showing a larger portion of the promoters downstream of the TSS.

The average profiles of each signal in Fig 2A are generally coherent, primarily illustrating the concentration of DNA accessibility in proximity to the TSS within the nucleosome-depleted region. However, to understand the local interaction between nucleosomes and TFs, it is crucial to examine the specific positions within the entire promoter region. Fig 2B shows examples of individual promoters and the relations between nucleosome binding, chromatin accessibility, and the potential nucleosome support described by the NF score.

The transcription factor NANOG is a general developmental factor involved in the proliferation and renewal of embryonic stem cells (46), and thus, it does not play a role in the differentiated cell lines. Fig 2B (left) displays the integrated profiles for the NANOG gene, which is not expected to be actively transcribed in either cell line. Therefore, the first significant observation is the lack of clear signals in both DNase-seq and H3K4me3 ChIP-seq data. In addition, the general accessibility in primary cell lines, described by the DHS score, is mainly zero overall, too, albeit the NF score is high near the TSS (highlighted region) with an NF score of 0.9–0.95. At the same location, a maximum of nucleosome occupancy (MNase-seq of 4.4 in K562 and 3.6 in GM12878) can be observed. The evolutionary conservation at this point is much higher with up to 2.9 compared with the near-zero level in the rest of the promoter. This indicates a regulatory nucleosome position that is currently blocked by a well-positioned nucleosome so that the activation of the NANOG gene requires an active displacement of this nucleosome.

In contrast, the transcription factor TP53 carries substantial importance in both the K562 and GM12878 cell lines, as it functions as a critical tumor suppressor (47). Thus, this gene stands as a representative case of a promoter subject to positive regulation. The profiles for this TF can be found in Fig 2B (right). The accessible regions are separated into two distant peaks on the level of DNase-seq and DHS score, one of which is located around the TSS and the other one at 800 bp upstream. Both accessible regions seem to be important in almost all primary cell lines (DHS score of 401 at the upstream peak and 398 downstream).

Both of these peak areas furthermore exhibit an absence of nucleosomes in these areas, evident in an MNase-seq signal of 0–1 at the upstream peak and of 0–0.3 in the downstream peak. Interestingly, the downstream peak is flanked by a high nucleosome occupancy of up to 1.5–2.1 to either side, potentially marking shifted nucleosomes and/or the well-described +1 and −1 nucleosome downstream of the TSS. The sequence support for nucleosome binding in the form of an NF score shows a broader plateau of 1 at the upstream portion of the accessible promoter and another maximum directly upstream of the TSS in the nucleosome-depleted region with a sharper profile and a short decline before it peaks again closely downstream of the TSS.

Fig 2A and B together leads us to say that locations in promoters where TFs compete with nucleosomes directly are the places with the highest nucleosome support. Whether these are occupied by nucleosomes is dependent on whether the gene is actually expressed. Thus, we are not of the opinion that all nucleosomes in the promoter are sequence-positioned, as suggested in the mean of all promoters in Fig 2A. In these overall promoter regions, the NF score stays constantly relatively high, but in the individual promoters in Fig 2B, it depends on the co-localization with competing TFs that mediate chromatin accessibility at local, defined places. The NF score exhibits a positive correlation with nucleosome occupancy in inactivated genes at regulatory key nucleosomes, where the nucleosome is deliberately positioned by sequence. In contrast, it correlates negatively once the corresponding nucleosome is displaced by targeted TF action at the same location.

Transcription factors

The regulation of gene expression is highly dependent on the binding of TFs. TFs exhibit a range of strategies in their competition with nucleosomes, displaying remarkable adaptability. Nevertheless, a common feature among them is that a region occupied by TFs should typically be devoid of nucleosomes. To give an overview of the nucleosome support for the diverse array of TFs, we calculated the mean NF score for all available ENCODE TFs at the peak positions of optimal IDR threshold peaks. The number of available files with at least two biological replicates is 151 for GM12878 and 391 for K562. Fig 3A shows a histogram of NF scores at each according individual peak point position as defined in the Methods section. It is notable that the overall mean nucleosome support of TFs is in general rather high with ∼0.71 and very consistent between the two cell lines with a difference within less than 1%. The difference in appearance between these two groups lies solely in the imbalance in the number of available experiments for each cell line. Furthermore, we calculated the mean NF score around the TF ChIP-seq peak location summarizing all of the TFs into one profile, which is depicted in Fig 3B. The peak of the NF score coincides precisely with the center of the experimentally validated ChIP-seq peaks. As expected, the NF score at this position matches precisely with the median NF score at the peak position depicted in Fig 3A’s histogram, measuring 0.71 at the peak’s maximum position. The figure illustrates that, in the vicinity of TF binding events, the MNase-seq signal experiences a consistent reduction (GM12878: signal decrease from 0.91 to 0.85; K562: signal decrease from 1.01 to 0.89). In contrast, the NF score reaches its peak precisely at these locations, suggesting the presence of sequence-intrinsic nucleosome positioning, particularly in regulatory regions that compete with initially well-positioned nucleosomes. The NF score lies in between what is to be expected for enhancer and promoter regions considering Fig 3. Because of the abundance of ChIP-seq peaks that are expected to be located in enhancer regions, the profiles might be biased toward enhancer TFs. To show whether promoters and enhancers behave differently, we added the information of H3K4me3 and H3K27ac signal for the TF ChIP-seq peaks in Fig S2. In this context, it is apparent that both signals exhibit a consistent depletion pattern in close proximity to the binding locations, demonstrating a strong similarity in their behavior.

Overall sequence-intrinsic nucleosome support around human TF binding locations.

(A) Sequence support at TF binding locations. Histogram of the NF score at the point-source peak position of all available ChIP-seq peaks on ENCODE for K562 and GM2878. The red line marks the mean of the NF score at peaks (GM12878: 0.7150, K562: 0.7114). (B) Average sequence-intrinsic nucleosome support (NF score) around the point-source peak of all TFs in comparison with nucleosome occupation (MNase-seq) in K562 and GM12878. It can clearly be seen that the sequence-based positioning signal is highest around the ChIP-seq peak location (0.71), whereas this is the position with the lowest nucleosome occupancy in vivo. This is evidence for the high importance of sequence to position nucleosomes based on the sequence at places where competition with TFs is highest.

TF profiles.

Extended figure relating to Fig 4B in the main article. Sequence-intrinsic nucleosome support, nucleosome occupancy, and histone modification signal for H3K4me3 and H3K27ac around the ChIP-seq peaks of all human TFs from ENCODE in hg19 as a mean profile surrounding the point-source peak. The sequence-intrinsic nucleosome support is highest around the TF binding locations. At the same time, the nucleosome occupancy signal is depleted. The profiles for both histone modifications are widely decreased around the ChIP-seq peaks, showing that both promoters and enhancers are expected to behave similarly in regard to the competition between TFs and nucleosomes. y-Axes: H3K27ac—ChIP-seq control normalized tag density; H3K4me3—ChIP-seq control normalized tag density; DNase—read depth normalized signal; MNase—density graph of signal enrichment; NF—NF score.

Because ChIP-seq only covers the final state of the displacement of the nucleosomes, we illustrate the relation of originally closed off chromatin and the sequence-based nucleosome positioning with a pioneer factor that serves the purpose of displacing the nucleosome and making the chromatin accessible. The pioneer transcription factor GATA3 has been described to bind well-positioned nucleosomes (48). To check whether this positioning can be seen on the DNA level as well, we analyzed 43,504 GATA3 binding sites with our global NF score and compared them with published nucleosome occupancy (MNase-seq) in MDA-MB-23 cells (49). Fig 4A shows the profiles of +500 bp around GATA3 binding sites, once in cells with GATA3 expression (+GATA3) and once with GATA3 knockout (−GATA3). It is evident that the area is supportive for nucleosome binding based on sequence as shown by an NF score that shows an increased peak of ∼0.15 around the binding sites from 0.35 to 0.5. In the case of cells without GATA3 expression, there is a shift of increased nucleosome occupancy (0.63–0.7) around the sites, showing increased nucleosome occupancy on top of the closed binding sites. However, upon GATA3 expression, the peak observed in the data undergoes a split into a double peak at 0.675, corresponding to both sites, while exhibiting a distinct valley at 0.62 directly at the TFBSs. This behavior suggests the presence of a robust nucleosome that was initially positioned based on the DNA sequence but subsequently displaced because of the influence of pioneer transcription factors.

Nucleosome formation (NF) score and nucleosome occupancy (MNase-seq) surrounding the transcription factors GATA3 and CTCF.

(A) GATA3. GATA3 binding sites of MDA-MB-23 cells with (green) and without (blue) GATA3 expression. The nucleosome occupancy is higher toward the binding sites in the cells without GATA3 expression. This peak is depleted in the GATA3 expression cells, and a double peak of nucleosome occupancy is left around the sites with a distance of ∼200 bp. The nucleosome support displays a pronounced peak precisely overlapping the TFBSs, indicating the presence of an intrinsically well-positioned nucleosome that is actively displaced by GATA3. (B) CTCF. CTCF ChIP-seq peaks from ENCODE in GM12878 and K562. The nucleosome occupancy shows a very distinct periodicity around CTCF peaks. This periodicity of nucleosome positions is not apparent on the level of sequence-intrinsic nucleosome support, which exhibits a single peak directly around the CTCF peak.

Insulators

Another category of regions with a rather high mean NF score of 0.66 in Fig 1A were insulator regions. Insulators are DNA regulatory elements that play a crucial role in organizing and maintaining the spatial organization of the genome. The protein CTCF is a key factor in mediating insulator function by facilitating chromatin looping and regulating gene expression (50). The nucleosome occupancy and sequence-intrinsic nucleosome support around 43,247 ChIP-seq peaks for CTCF in K562 and 43,631 peaks in GM12878 are compared in detail in Fig 4B.

The nucleosome occupancy shows a very stable, phased array of nucleosomes around the CTCF peaks with an ∼200-bp periodicity and an ongoing MNase-seq signal decrease per peak from 1.8 to 1.2 with increasing distance to the TSS. Although the nucleosome occupancy is lowest directly at the CTCF peaks, the sequence-intrinsic nucleosome support is rather high specifically in these regions with an NF score of about 0.77 at the peak. The observed periodicity in nucleosome occupancy is completely absent on the sequence level, indicating a strong influence on nucleosome positioning at CTCF sites but active energy-driven chromatin remodeling once the site was cleared and the factor is bound to the DNA.

Sequence-intrinsic nucleosome support in transcriptional regulation

In the previous part, we focused mainly on the role of DNA for nucleosome positioning in promoters, in proximity to the TSS and in competition with TFs. The following sections will explore other functional elements of the genome along the transcriptional axis in a similar way.

Active chromatin remodeling upon transcription initiation

It is known that there are a multitude of influences on the dynamic positioning of nucleosomes. One such phenomenon is the phasing of ordered nucleosome arrays downstream of the TSS upon transcription initiation (51). Fig 5A shows the comparison of nucleosome occupancy (MNase-seq) and the DNA-intrinsic nucleosome support (NF) around the TSS of all human RefSeq genes. The proposed periodic downstream array of nucleosome positions is apparent, marked by arrows, stretching with a periodicity of 200 bp at MNase-seq levels between 0.75 and 0.9 (in GM12878) or 0.65 and 0.85 (in K562), respectively. At the same time, although the NF score is generally still relatively high close to the TSS and the +1 nucleosome with a value of 0.8–0.85, the periodicity is not mirrored on the sequence level and the general sequence-intrinsic nucleosome support decreases further downstream. Therefore, it is reasonable to conclude that the dynamic relocalization of nucleosomes because of ATP-dependent remodeling can overwrite the sequence-mediated nucleosome guidance.

MNase-seq and NF score for different elements from transcription initiation to transcription termination.

(A) Promoters. The profiles are consistently aligned to ensure that the promoter located at the 5′ end is consistently positioned on the left side. The central vertical line represents the TSS location. The nucleosome occupancy declines around the TSS and shows an apparent periodicity in both directions. The NF score rises in anti-correlation toward the TSS and does not indicate the proposed periodicity. This implies there is support for nucleosome binding for the NDR, which is accessed by TFs but no sequence support that places the periodically following nucleosomes in the downstream array. (B) Exons. The left image is centered specifically on the intron/exon (IE) junction, whereas the image on the right aligns the exons precisely over their midpoints. The NF score shows high sequence-intrinsic nucleosome support and an experimentally high nucleosome occupancy directly downstream of the exon starts, which can be observed even more stably over exon centers. This correlation strongly suggests the functional significance of the observed sequence-defined nucleosome. (C) Poly(A). There is a decline in both sequence-intrinsic nucleosome support (NF score) and nucleosome occupancy (MNase-seq) around the transcription termination sites. This behavior indicates nucleosome-repelling positions, encoded in DNA sequence.

Sequence-intrinsic nucleosome occupation potential of exon regions

Although it was just shown that there is a sequence-independent creation of nucleosome arrays downstream of the TSS, that does not entirely describe the sequence-intrinsic local nucleosome support for other regions along transcriptional regulation. It is described in literature that exons are targets of a particularly high nucleosome occupancy, which is related to transcription speed of RNAPII and the recognition of alternative splicing (11). To test whether this principle is also evident on the nucleotide sequence level, we analyzed the exons of all human genes by comparing the MNase-seq signal for K562 and GM12878 and the NF score again. Fig 5B shows both average signals centered around 267,864 exons with strand orientation taken into account. The figure on the left employs an alignment centered around the intron/exon (IE) junction, specifically focusing on the exon starting points. Conversely, the right-hand figure emphasizes the alignment based on the central exon midpoints. Generally, the figures indicate that there is a nucleosome positioned over the center of the exon. The midpoint-centered figure shows that the nucleosome occupancy and the NF score show a distinct peak over the exon regions. The extent of the rise of the NF score toward the exon center is as high as ∼20% from slightly above 50% to almost 70% and is distributed mainly between 100 bp upstream and 150 downstream of the exon midpoints. The nucleosome occupancy mirrors this pattern of a central peak of about 0.3–0.4 from 0.9 up to 1.25 in both cell lines around the exon. In addition, a notable observation indicates a relative depletion of MNase-seq signal directly upstream of the peak, which is more clear when aligning the exons over the IE junction. In this view, it is also apparent that the main nucleosome is positioned directly downstream of the IE junction. That the middle-centered nucleosome can be understood as equivalent to the nucleosome directly downstream of the exon start can be explained with the median of the RefSeq exon length being 129, thus being not large enough to fit multiple nucleosomes. The general tendency of a single sequence-positioned nucleosome in the center of an exon holds true when separating the exons by length. For a more detailed analysis of these patterns, Figs S3 and S4 provide an enhanced representation by separating the exons based on their length, again around the midpoints and around the IE junction, respectively. The images shown here encompass all human exons. Because the first exon is essentially the TSS, Fig S5 shows the same profiles, once for all exons except the first of each gene and once only for the first exons. Because of their comparatively small amount, the principal image does not change when filtering out the first exons. The isolated first exons, however, do not exhibit the same occupation of exons but rather show the same nucleosome positioning as shown in the TSS analysis in Fig 5A with the clear +1 nucleosome.

Nucleosomes around exon midpoints.

MNase-seq and NF scores for all human exons, centered around the center of all exons. The quantiles separate different lengths of exons. The lower quantiles denominate shorter exons. In particular, they denominate these exon lengths: 1 = 1–90 bp (median 69 bp), 2 = 90–129 bp (median 109 bp), 3 = 129–190 bp (median 153 bp), 4 = 190–91,671 bp (median 333 bp). The strand orientation is taken into account and always 5′-3′–directed. The NF score shows a high sequence-defined nucleosome support for exon regions, and experimentally, high nucleosome occupancy can be observed. The longer exons generally display a higher NF score.

Nucleosomes around exon starts.

MNase-seq and NF scores for all human exons, centered around the intron/exon junctions. The quantiles separate different lengths of exons. The lower quantiles denominate shorter exons. In particular, they denominate these exon lengths: 1 = 1–90 bp (median 69 bp), 2 = 90–129 bp (median 109 bp), 3 = 129–190 bp (median 153 bp), 4 = 190–91,671 bp (median 333 bp). The strand orientation is taken into account and always 5′-3′–directed. The nucleosome occupancy is highest downstream of the TSS for all exon lengths. The IE junction is based on the apex of the rise of nucleosome positioning from an upstream depletion, around 75 bp in distance, toward a clearly positioned nucleosome 75 bp downstream.

First exons.

MNase-seq and NF scores for all human exons except the first (top) and for only the first exons, essentially the TSS (bottom). It is evident that the exons excluding all first exons are characterized by nucleosome occupancy over their midpoints, whereas first exons show only a clearly defined downstream nucleosome, which can be interpreted as the +1 nucleosome.

Sequence-dependent nucleosome occupancy prevention around poly(A) sites in transcription termination

Another important part of gene regulation and splicing is the nucleosome occupancy at poly(A) sites as the place for transcription termination. These were shown to be explicitly free of nucleosomes, thus making them accessible for proteins marking the sites and interacting with RNAPII in transcriptional termination (13). We obtained the ENCODE MNase-seq signal for K562 and GM12878, as well as the NF score, for a list of 570,740 poly(A) sites in the human genome (52). These are classified by the original authors into clusters of biological origin. All sites are oriented by strand in 5′-3′ direction.

Fig 5C shows the nucleosome occupancy and NF score for all poly(A) sites in the mapping together. There appears to be a direct positive relation between the sequence-intrinsic nucleosome support and the nucleosome occupancy. Both nucleosome occupancy and sequence-intrinsic nucleosome support exhibit a decrease around the sites in the central 250 bp. There is a 20% drop in nucleosome support relative to the already low level of 35%. It can be noted that the sequence support for nucleosome binding seems to be slightly higher in the upstream direction. The same overall decrease is apparent on nucleosome occupancy in both cell lines with a drop of 0.25–0.3 in MNase-seq signal. A proposed increase in nucleosome occupancy downstream of the poly(A) site in vivo relative to the upstream direction (13) cannot be observed overall. Nevertheless, outcomes for specific subgroups of poly(A) sites can differ. Fig S6 gives an overview over all poly(A) sites separated by clusters defined in the PolyASite 2.0 atlas. There, it can be observed that the general repelling trend can be observed for all subclusters, except for sites that coincide with exon regions, which exhibit a similar pattern as shown above for exons, that is, a sequence-positioned nucleosome over the exon. The aforementioned difference of a higher nucleosome occupancy downstream of the site can be shown there for the DS cluster, which denotes the TTS 1 kb downstream of a terminal exon. In contrast, the higher NF score in the mean profile upstream of the sites can most probably be explained by the group of sites located in the AU group, in which the signal is higher because of its proximity to the TSS, which has been demonstrated to have a rather high NF score in Fig 5A.

Poly(A) sites.

Detailed overview over all poly(A) sites from the PolyASite 2.0 atlas, separated by cluster (AE = anti-sense to an exon; AI = anti-sense to an intron; AU = 1,000 nt upstream in anti-sense direction of a transcription start site; DS = 1,000 nt downstream of an annotated terminal exon; EX = exonic; IG = intergenic; IN = intronic; TE = terminal exons) (52). All groups apart from exonic poly(A) sites exhibit a local decrease around the site in both sequence-defined nucleosome support (NF) and nucleosome occupancy in both cell lines K562 and GM12878.

Low-frequency differences in HelT are the most discriminative sequence feature for nucleosome positioning

The random forest was selected for modeling the classifier underlying the NF score because of its intrinsic ability to generate feature importance scores. The input features consist of the power spectral density (PSD) of frequencies inherent in underlying DNAshape structures (see the Classifier section in the Materials and Methods section). In Table 1, the top 10 input features for distinguishing between nucleosomal and linker DNA are listed. The mean decrease in Gini impurity is given for all frequencies of the PSD-transformed DNAshape values. It is evident that most of the best input features to distinguish between nucleosomal and linker DNA are low-frequency differences in HelT. In the original DNA structure space, low frequencies indicate long-range structural changes rather than short local perturbations.

Table 1.

Feature importance of the top input features for the random forest classifier as a mean decrease in Gini impurity.

As HelT denotes the rotational angle between two adjacent bp along the helix direction, this effect can be understood as what happens when taking both ends of a 147-bp DNA double strand and twisting them in opposite directions. Fig S7 shows the visualization of the mean profiles for each DNAshape feature, as well as the PSD-transformed information for both nucleosomal and linker DNA from the training set, and underlines the differences between both groups in the low-frequency structure changes visually. It has been shown in this study that the NF score is highest at regulatory meaningful positions, which are targeted by TFs. It is worth noting that the long-range structural disparities between nucleosomal and linker DNA contrast with the localized and short DNA structure features known to be recognized by TFs at their binding sites (53). Therefore, attraction for the binding of nucleosomes and for TFs is encoded at the same locations but on different frequency ranges.

Visualization of the training data.

Left: average profiles of DNAshape values over all training samples from the Guo et al set of linker versus nucleosomal sequences. Right: average profiles of the power spectral density for each of the DNAshape features. The datasets appear separable in principle. The power spectral density profiles are more different between both groups on low frequencies.

Discussion

In this study, we analyzed the role of DNA sequence for the positioning of nucleosomes. For that purpose, we developed a random forest classifier to distinguish nucleosomal and linker DNA and used it to derive a score to estimate the sequence-intrinsic nucleosome support of different functional elements in the human genome. The degree to which nucleosome positioning relies on DNA sequence has been the subject of extensive debate (16). These questions cannot be answered by analyzing either chromatin accessibility experiments or DNA sequence individually. On the one hand, the experimental verification of chromatin accessibility for any given cell type will only describe the result of the multitude of counteracting influences on nucleosome positions. On the other hand, DNA can explain the importance of nucleotide sequence on their potentially favored binding locations, without reliably predicting the actual nucleosome occupancy in any given cell type in vivo. Therefore, we compared the results of our sequence-intrinsic nucleosome support in the form of the NF score with experimental data for chromatin accessibility and histone modifications and described the possible combinations of in vivo nucleosome positioning with its underlying DNA sequence support. We could demonstrate that DNA-intrinsic nucleosome support is not directly equivalent to the resulting nucleosome occupancy in a specific cell line in vivo. Rather, there is a positive influence on nucleosome binding in the DNA sequence of genomic regions, which are involved in actively regulating transcriptional processes, such as promoters, enhancers, exons, and in the center of insulators. Simultaneously, certain locations exhibit lower nucleosome support, which can serve functional purposes, such as well-positioned nucleosome occupancy arrays downstream of the TSS and in the vicinity of insulators. Conversely, there are instances where lower nucleosome support indicates non-functionality, as observed in constitutive heterochromatin. That sequence-positioned nucleosomes in regulatory regions are a functional effect and not the result of experimental biases in relation to linker sequences in the underlying training data from Guo et al is illustrated by Fig S8. This image compares the training sequences in terms of their relative proportion over different biological contexts. There are no major differences in the origin of both categories, which are mainly distributed over the most abundant regions of the human genome such as distal intergenic regions and introns.

Training sequence distribution over biological contexts.

Distribution of Guo et al (2014) training set samples in biological settings. It is evident that both nucleosomal and linker sequences are found predominantly in the abundant distal intergenic and intronic locations. There is no experimental bias that would explain the principal difference between both categories to cover regulatory elements (nucleosomal) or randomly positioned nucleosomes (linker). Principally, there are slightly more nucleosomal sequences found in promoter regions than in the linker group. However, given the relatively low occurrence in both groups, this is within the range to be expected under our hypothesis of a meaningful enrichment of nucleosomal sequences within regulatory important locations. To obtain the distribution, we used the ChIPseeker R package (54) to categorize the Supplementary BED files SD1 and SD2 that contain the hg19 mapping of the original FASTA files from Guo et al (2014).

The increased reliance on DNA sequence in gene regulatory regions could be shown in this study from a variety of biological angles. The apparent general trend in promoter regions could be narrowed down to be due to the direct places of competition with TFs. The nucleosome support by sequence could be attributed to evolutionary pressures favoring the development of DNA sequences that facilitate optimal nucleosome binding, effectively regulating access to these regions. This effect has been described to be particularly important for organisms of higher complexity (55). Especially, these highly sequence-intrinsic and evolutionary conserved positions can make a large difference in the regulation of gene activity, and thus, the sequence can help to increase the resistance to accidentally triggering gene expression. The counteracting process to this state is an active eviction of the blocking nucleosome by a combination of (pioneer) TFs, histone modification enzymes, and active chromatin remodeling to gain accessibility for a particular position (56). Because of the inclusion of the relevant TFBSs within the confined nucleosome positions, it can be inferred that there are a presumed co-evolution of proficient nucleosome organization and the existence of high-quality binding sites for activating TFs. Therefore, we interpret the crucial function of DNA in these “regulatory” nucleosomes as a mechanism to regulate the competition between nucleosomes and TFs for spatial availability. This competition affects the energy required to activate biological processes, allowing for dynamic control and modulation. This co-localized competition raises the question of how sequence-encoded nucleosome preferences and well-defined TFBSs can evolve at the same positions. The answer lies in the analysis of feature importance derived from the random forest model. This study demonstrates that the primary predictors for distinguishing between nucleosomal and linker DNA stem from low-frequency, that is, long-range, differences in HelT. Simultaneously, it is well documented that TFs recognize very local, high-resolution, that is, high-frequency, DNA structure signals (53). This leads to the assumption that the nucleosome attraction and the counteracting TFs are encoded at the same positions but on different frequencies of DNA structure changes. Further evidence for the systematic co-localization of regulatory proteins (high DHS score) with highly nucleosome-attracting regions (high NF score) can be observed in Fig S9. That the accessibility of these sequence-defined nucleosomes plays a significant role in regulating gene expression is further supported by Fig S10, in which gene expression values are compared between groups of genes based on a correlation or anti-correlation of the NF score and chromatin accessibility data. The role of DNA sequence for different kinds of TFs under this assumption has been partially analyzed for different TF families (57, 58). Nevertheless, it would be worthwhile to research the differences in sequence-intrinsic nucleosome support in regard to pioneer functionality. The co-evolution to occupy highly important places would suggest a high nucleosome support at pioneer TFBSs. However, this presents a notable challenge as the binding sites of pioneer transcription factors are often located in close proximity to the subsequent TFs that bind in a consequent step and that are not directly in competition with the nucleosome. Furthermore, it is well established that pioneer TFs can employ various modes of nucleosome eviction, which may not occur directly at the target site of chromatin accessibility (7, 59). Finally, although there is a difference in the set of TFs that are responsible for regulating the chromatin accessibility of promoters and enhancers (60), we could show that for the general relationship of nucleosome positioning in the competition with TF binding, both groups show a similar behavior.

Sequence-defined nucleosome support in accessible genomic regions.

To show the relation between highly important gene regulatory regions and the sequence-intrinsic nucleosome support, we compared the NF score and the DHS score (see the Materials and Methods section) with the Fisher test. For that purpose, we took the individual base pairs of all human promoters and divided both measures into two bins each, according to their median. The median for the high-resolution NF score is at 0.76 and the DHS score at 1 (i.e., 1 of 403 primary cell lines is open at that particular bp). There is a significant (Fisher test, P-value < 2.2 × 10⁻¹⁶, odds ratio: 4.416, 95% confidence interval: [4.409 4.423]) relation between regions with high sequence-intrinsic nucleosome support and accessible chromatin in vivo, indicating a co-evolution of nucleosome and transcription factor attraction in human promoters.

Sequence-positioned nucleosomes affect gene expression.

We calculated the Spearman correlation between the high-resolution NF score in all human promoters with the respective MNase-seq and DNase-seq signal per promoter. Then, we tested whether the gene expression is lower when the NF score and MNase-seq showed a significant positive correlation versus a negative correlation. Respectively, we tested for a higher gene expression in genes with a positive correlation of the NF score and DNase-seq versus gene expression in negatively correlated genes. The table shows the P-values of the Wilcoxon test per cell line (K562, GM12878) and comparison (MNase-seq, test direction “less”; DNase-seq, test direction “greater”). The P-values are adjusted with the Bonferroni–Holm method. The significant P-values suggest that there is an influence on gene expression based on the displacement of sequence-intrinsically positioned, “regulatory” nucleosomes. The definition of promoters follows the methods of the main article. The gene expression values are obtained from ENCODE (GM12878: ENCFF873VWU [rep1], ENCFF345SHY [rep2]; K562: ENCFF928NYA [rep1], ENCFF003XKT [rep2]) (41)

In contrast to these nucleosome positions determined by DNA sequence, constitutive heterochromatin, although commonly occupied by nucleosomes, does not demonstrate a comparable reliance on DNA sequence for precise nucleosome positioning. In analogy to the previous explanation of DNA evolution toward a targeted occlusion of TFBSs to prevent random gene activation, we propose a lower relevance of sequence for heterochromatin. The higher order packing of dense chromatin is a more unspecific process in which the accessibility of small individual positions on the genome is not expected to have such severe consequences as the revealing of TFBSs in cis-regulatory elements.

Apart from the evolution of nucleosome positioning advantage at highly competitive sites, we found other instances of locations that appear to have high sequence-intrinsic nucleosome support. In particular, exons show considerable nucleosome support relative to their surroundings. It is well described that the occupation of exon regions with nucleosomes is important in the pausing of RNAPII to ensure the correct transcription of a gene, as well as the marking and spatial relocation of exons for the splicing machinery (11, 61). Certain studies propose that the sequence specificity for nucleosomes at exons is primarily encoded at the exon’s outset, specifically in the form of a specific dinucleotide distribution at the intron/exon junction (62). However, the positioning of the nucleosome based on the sequence is found in this study directly downstream of the exon start. This finding is in line with the theory that the nucleosome in the exon stalls the RNAPII just upstream of the IE junction (63). The sequence-encoded attraction of nucleosomes to exon locations might therefore be another example of reducing the energy needed to guide nucleosomes to consistently regulate RNAPII dynamics.

The principle of using the DNA to encode a favorable structure to guide nucleosomes to exon regions might again be an efficient means of increasing the intrinsic nucleosome support. In contrast to strengthening the nucleosome in competition with TFs, the nucleosome at exons might simply be guided to its place of effect without the need for further displacement by other processes. The difference in the non–sequence-dependent nucleosome positioning at the TSS could be grounded in the reduced energy needed to stabilize nucleosomes permanently at exons compared with the potential for dynamic shifting of nucleosomes around TSS upon gene activation. This provides another striking example of putative co-evolution that is directed to attract nucleosomes on the sequence, which underlies a drive to reliably encode proteins at the same time.

The use of sequence to influence nucleosome positioning is evident not only in its ability to positively guide nucleosomes to specific locations. We showed the opposite effect for poly(A) sites, marking the end of transcription. The substantial contrast in sequence-intrinsic nucleosome support between poly(A) sites and the surrounding DNA clearly indicates the strategic prevention of nucleosome binding at these specific locations. The reason for this feature could be the attempt to keep chromatin at the TTS constantly accessible to easily allow the binding of proteins that form complexes together with the RNAPII to terminate transcription, cleave the mRNA, and do the polyadenylation in the consecutive step. The nucleosome rejection fits well with the general observation of low nucleosome occupancy over poly(A) sites. However, there are hints about the role of increased nucleosome occupancy around one end of the sites, regulating alternative polyadenylation and depending on the status of active transcription (13). This could be the reason for a difference in nucleosome support toward both directions surrounding poly(A) sites in terminal exons. Here, we observed slightly lower nucleosome support downstream of poly(A) sites for terminal exons but not as a general trend over all groups. Also, we could show an increased nucleosome support upstream of the TSS for the group of sites downstream of the TSS, in which the increase is easily explainable by the overlap with the generally sequence-defined promoter region. It cannot definitively be determined whether the otherwise relatively constant nucleosome occupancy up- and downstream arises from technical aspects in the normalization of ENCODE MNase-seq data or whether it is attributed to variations in the splicing patterns of different clusters of poly(A) sites, particularly when situated in regions other than terminal exons, or whether all of these sites possess transcription termination functionality at all.

Opposed to the steering of functional nucleosomes by nucleotide sequence, we identified processes that exhibit strong nucleosome positioning that does not manifest on the DNA sequence level. Notably, downstream of the TSS, a well-established pattern of regularly spaced nucleosome array formation is observed. At these locations, the nucleosome support constantly decreases toward intronic locations without mirroring the periodicity of the nucleosome occupancy. This process is known to be mediated by active chromatin remodeling complexes, which are organizing nucleosomes into an evenly spaced array using energy from ATP (10). This process serves as an example of an important transcriptional mechanism, which is not directly supported by local DNA structure. Regarding the nucleosome arrays adjacent to the TSS, it can be speculated that the +1 and −1 nucleosomes are influenced by sequence-specific influences. This assertion is supported by the relatively high NF score, aligning with studies showing the special significance of these particular nucleosomes (8, 64). These studies suggest a discernible dinucleotide enrichment in proximity to these nucleosome locations, implying a mechanism of “statistical positioning” for the co-occurring nucleosome arrays. However, the proximity of these nucleosomes to the main areas of TF–nucleosome competition in the promoter could easily lead to misinterpretation of their reliance on DNA sequence alone.

In conclusion, it is worth highlighting that these processes of sequence-intrinsic and sequence-independent nucleosome positioning may not necessarily be mutually exclusive. This phenomenon is particularly pronounced at insulator sites, where the transcription factor CTCF plays a crucial role in delineating TAD boundaries. Here, it could be observed that on top of the actual binding location there is very strong sequence-defined support for the nucleosome, suggesting the previously described mechanism of competition and co-evolution of TF–nucleosome interaction. Yet, as a result in the surroundings of CTCF-bound locations, there is a clear well-positioned nucleosome array, which is the result of active chromatin remodeling without the help of sequence support. Although the binding of CTCF underlies the aforementioned competition with positioned nucleosomes, the spacing and symmetry of the surrounding array are rather linked to active chromatin remodelers, instead (5).

This is just one example of the very close, local interplay of sequence support and energy-dependent remodeling of chromatin to decide the final in vivo positions of nucleosomes in transcriptional regulation. It is impossible to favor any of these principles in their importance for the collaborative regulation of transcription, from the formation of TAD boundaries, over the regulation of gene activity through the accessibility of chromatin for TFs or support mechanisms of the tasks of RNAPII, to the clean termination of transcription, thus marking the last step in which the DNA is involved directly.

Materials and Methods

Datasets

The experimental chromatin accessibility datasets were obtained from ENCODE for the human genome in hg19 (41). They comprise MNase-seq data as a density graph of signal enrichment, described in reference 65; DNase-seq data as read depth normalized signal, described in reference 66; and peak-called histone modification ChIP-seq data as bed narrowPeak locations for K562 and GM12878 and raw signal for H3K4me3 and H3K27ac as control normalized tag density, described in reference 67. Besides, we used CTCF ChIP-seq peaks from ENCODE for our insulator analysis. Furthermore, we used the comprehensive list of ChIP-seq narrowPeak locations for all TFs that are available on ENCODE for hg19. When referring to the central peak position related to the ENCODE data, we reference the point-source peak defined in the narrowPeak format standard. Because of their extensive length, we show the ENCODE IDs for all the aforementioned datasets in Tables S2, S3, and S4.

Table S2 ENCODE IDs for individually mentioned datasets used in this study.

Table S3 ENCODE IDs of all ChIP-seq experiments used for the TF analysis in GM12878.

Table S4 ENCODE IDs of all ChIP-seq experiments used for the TF analysis in K562.

In addition, we used nucleosome occupancy data to evaluate the influence of the GATA3 TF in MDA-MB-23 cells with and without GATA3 expression. The according MNase-seq data were taken from Takaku et al (2016) (49). We calculated the mean of all three replicates for the condition of GATA3 expression and for cells without GATA3 expression.

Furthermore, we evaluated the evolutionary conservation at a base pair resolution with the PhyloP score (68), which was obtained using the UCSC Table Browser (45).

We compared all these data along a variety of different genome locations. The broad overview over different biological contexts is defined by chromHMM for K562 and GM12878 and defines 15 meaningful states originally derived from histone modification ChIP-seq data (44).

The gene mapping we use is based on RefSeq in hg19 (69) and was also obtained with the UCSC Table Browser. We exclude genes from chromosome Y, because not all experimental chromatin accessibility data are available for it. In particular, we use the remaining 30,142 unique TSS locations to define genes and promoters. Our promoter definition is defined as within 1-kb upstream of each TSS defined by RefSeq. In case of multiple alternative TSS locations, all of these were considered as individual gene starts. Besides using RefSeq for TSS, we also extracted all exon locations from the same mapping.

Moreover, we analyzed the elements at TTSs in the form of poly(A) sites. The poly(A) site mapping with a separation into clusters of biological origin is taken from PolyASite 2.0 (52). To obtain single locations for the center of each site, we used the representative cluster location from the unique cluster ID column. The original assembly of the data is hg38. We used liftOver to convert the coordinates to hg19 via the rtracklayer R package (70).

Software

Apart from the individual software that is referenced in its due place, the analyses were performed using R (v. 4.3) (71), and we used the packages rtracklayer (v. 1.60.0) (70), ggplot2 (v. 3.4.2) (72), GenomicRanges (v. 1.52.0) (73), ROCR (v.1.0.11) (74), and randomForest (v. 4.7.1.1) (75).

Classifier

The classifier is based on the human Guo et al (2014) dataset (33), which is in turn derived from an in vivo mapping of nucleosome positions in CD4⁺ T cells (39). It contains 4,573 non-redundant FASTA sequences, which are categorized as either nucleosomal (n = 2,273) or linker sequences (n = 2,300). A mapping of these sequences to the hg19 human reference genome is in addition provided as BED files (see the Data Availability section). The labeling strategy of the original authors is based on a ranking of sequences by a summation of the underlying nucleosome occupancy counts. Therefore, nucleosomal sequences are the most consistently occupied genome regions, whereas linker sequences are understood as stretches of DNA that do not show a nucleosome occupancy signal. Because nucleosomes usually bind 147 bp of DNA, the same length is used for the linker sequences, although depending on the definition, these can show different lengths in the genome. Fig 6 shows the workflow for setting up the classifier (Fig 6A), transforming the input data and its genome-wide application with a sliding-window approach (Fig 6B). The FASTA files are transformed with the DNAshapeR tool (v. 1.28.0) (43), resulting in a set of numerical vectors for Roll, HelT, ProT, and MGW with one value per bp each. The mean profiles of each DNAshape feature are compared between nucleosomal and linker sequences in Fig S7 (left). In addition, a power spectrum (PSD) is applied to the DNAshape values in order to strengthen positional independence of the periodic nucleotide patterns, emphasizing the assumed underlying periodic sequence pattern and making it suitable for the use in a random forest classifier. That step uses the spectrum function in R with a taper of 0.1 and a span of 10. In the next step, the resulting spectral densities are concatenated into a feature matrix, which is the input for a random forest classifier. The input features for the classifier consist of the PSD values per frequency for each DNAshape feature. Inside the model, features are denoted by a combination of the DNAshape feature and frequency. For instance, the input feature HelT_5 represents the power of the underlying HelT sequence at a frequency of 5 per 1/147 bp. A high power for low-frequency numbers indicates a significant contribution from long-range signals. Fig S7 visualizes the PSD profiles for nucleosomal and linker DNA. The random forest classifier was trained with the randomForest package in R with default parameters. The result of the classification is the categorization of a sequence being either nucleosomal or linker and serves as the basis for the two types of NF scores that will be described further in the NF score section. The particular final number of input features at the end of the transformation is 288. This results from the basic sequence lengths (147 bp), the number of DNAshape parameters used for training (4), and a subtraction of some missing DNAshape values at the edges of the sequence where the centered k-mer mask is incomplete (12). These 576 features are halved because of the redundant nature of both sides of a full PSD spectrum to the final number of 288 input features.

Processing and classification of DNA sequences—deriving the NF score.

(A) Random forest model is used to annotate DNA sequences as nucleosomal or linker. To build the feature matrix, DNAshape features are obtained for the input sequences and further transformed into a frequency power spectrum. These PSD profiles are then concatenated into the final input feature matrix to use in the random forest model for training and classification. (B, C) To create a pre-calculated mapping for the human genome, we applied the classifier genome-wide in a sliding-window approach with a step size of 50 bp and a detailed version for promoters with a step size of 7 bp. Both resources are provided as a downloadable bigWig file (see the Data Availability section).

The performance of the random forest classifier based on 10-fold cross-validation is shown in Table 2. The accuracy is measured with 82.01% (±1.55) and the AUC with 0.89 (±0.01). Furthermore, we compiled a range of different random forest models based on all of the previously mentioned data pre-processing steps. Table S5 lists performance measures for classifiers trained on raw sequences, pure DNAshape values, and the final PSD frequency transformation. It demonstrates that with each further transformation, there is a substantial gain in performance, which underlines the suitability of the frequency modeling to cover the underlying periodicity proposed by Segal et al (2006) (20).

Table 2.

Performance measures for the random forest classifier, evaluated with 10-fold cross-validation.

Table S5 Performance for different stages of training data pre-processing, evaluated with 10-fold cross-validation. The numbers in parentheses denote the SD of these measures within the 10-fold. Each transformation step exhibits a gain in performance.

The most recent classifiers use CNNs and reach accuracies of up to 0.889 (CORENup (76)) and an AUC of up to 0.94 (LeNup (35)) on the given dataset. However, although CNNs are undoubtedly valuable for such prediction tasks in principle, we opted not to use them on this specific dataset because of their propensity for overfitting on smaller datasets. Furthermore, the random forest methodology enables us intrinsically to derive feature importance scores from the model and we have a direct connection between the NF score resource and what determines the score on the basis of periodic DNA structure changes in the form of PSD. The feature importance was obtained using the importance function from the randomForest package in R and is provided as a mean decrease in Gini impurity. The ranking of each input feature by importance can be viewed in Table 1.

The present in vivo training dataset is a common standard in the prediction of nucleosome positioning. We want to point out that the use of in vivo data has the potential to introduce bias in predicting nucleosome positioning events, independent of the chosen classification method. Although previous studies have demonstrated the promise of incorporating in vitro nucleosome maps (55, 77), we find it noteworthy that despite training on in vivo data, we can show that the classifier is identifying the underlying sequence importance rather than the in vivo nucleosome occupancy that is resulting from the dynamic displacement of sequence-positioned nucleosomes in living cells. This is especially evident in the first section of the results about the counteracting relationship of high sequence-intrinsic nucleosome support in regulatory important regions against the low actual nucleosome occupancy in vivo at these positions of nucleosome–TF competition.

Score calculation

Nucleosome formation (NF) score

The classifier described above gives a binary prediction for a DNA segment of 147 bp in length. To allow for easier handling of the estimates of nucleosome support based on the nucleotide sequence, we pre-calculated two sets of scores for the human genome in hg19 assembly with different local resolutions as seen in Fig 6C. (1) The low-resolution score is a simple application of the classifier in a sliding window to all regions of the human genome with a fragment input length of 147 bp and a step size of 50 bp. The resulting NF score is an estimate of either 1 (nucleosomal) or 0 (linker) for the central 50 bp in any given 147-bp bin. (2) The high-resolution NF score provides a continuous score between 0 and 1 and is based on a much more local sliding window with a segment size of 147 bp and a step size of 7 bp. All base pairs are thus part of multiple (147/7 = 21) predictions, and the score represents the relative frequency of nucleosomal prediction results that this particular bp has been part of. To give an example, if an individual bp has been part of 21 sliding windows, of which 14 were classified as nucleosomal and 7 as linker, then the score at that base pair is 14/21 = 0.67. Because of the introduction of artifacts, otherwise the step size needs to be a divisor of the window size of 147 bp. Thus, we chose a resolution of 7 bp as it is the most local one next to a resolution of 1 bp.

Because of resource efficiency, the high-resolution score was only pre-calculated for all human promoter regions, because we want to give a good local estimate of the importance of the sequence for nucleosome positioning in specific individual promoter regions. In contrast to the high-resolution score, the step size of 50 for the low-resolution score is purely pragmatic because there is no further integration of the classification results. If not specified otherwise, we always use the low-resolution score, because we normally show the mean of nucleosome support profiles over a large number of fragments and there is virtually no difference between averaging the low-resolution 0/1 predictions and doing the same with the 0–1 high-resolution score. To illustrate the sufficiency of the 50-bp low-resolution score when examining larger sets of data in an overview fashion, we sampled 10,000 of the promoters used in this study randomly and calculated mean profiles with both the high- and the low-resolution scores and compared their similarity in Fig S11. The difference between both profiles varies maximally between 0.01 and 0.015. Both resources are easily accessible as bigWig files and listed in the Data Availability section. In addition, the code to classify nucleosomal versus linker DNA can be found there. The NF score described in this section is what is in this work referred to as sequence-intrinsic nucleosome support.

NF score resolution comparison.

We calculated and provided two resources. The high-resolution NF score is calculated for all human promoters and is built by calculating the relative frequency of nucleosomal DNA predictions with the classifier for overlapping sliding windows with a step size of 7 with the result of a score between 0 and 1. In contrast, the low-resolution score is created by applying the classifier in a binary way so that each window yields an array of 50 central bp within a 147-bp window of either 0 (linker) or 1 (nucleosomal). The profiles show the mean score distribution of 10,000 randomly sampled promoters from RefSeq. The vertical line at 0 bp denotes the TSS.

DHS score

In this study, we analyze the influence of sequence on nucleosome positioning, However, the resulting accessibility in particular cell lines can be very cell type–specific. Therefore, we created a score to generalize the chromatin accessibility over multiple cell lines. The DHS (DNase hypersensitive sites) score is a measure to estimate the generalized accessibility of any given genome location. For that purpose, we obtained peak-called DNase-seq data for all 403 primary cell lines from ENCODE in humans. These were integrated as a simple count of how many primary cell lines show a DNase-seq peak at any given base pair.

To give an example, a score of 403 at a specific bp means that all of the available cell lines show a DNase-seq peak overlapping that position. A score of 5 means that only five specific cell lines show accessible chromatin in the form of a DNase-seq peak at the specified bp. In the study, the score is mostly used in the context of indicating TF binding.

Data Availability

The code to execute predictions can be found at https://gitlab.gwdg.de/MedBioinf/generegulation/nfclassifier. The NF scores in bigWig file format can be downloaded at https://bioinformatics.umg.eu/resources/nfscore. The Guo et al training dataset was further mapped to the human reference genome hg19 in the form of BED files using the Rsubread R package (78). These mappings are appended as additional Supplemental Data 1 and 2, separated by nucleosomal and linker sequences.

Supplemental Data 1.

BED Mapping of nucleosomal training sequences (hg19).

Supplemental Data 2.

BED Mapping of linker training sequences (hg19).

Acknowledgements

We thank Argyris Papantonis and Alexander Ecker for fruitful discussions. M Sahrhage is a member of the International Max Planck Research School for Genome Science (IMPRS-GS), part of the Göttingen Graduate Center for Neurosciences, Biophysics, and Molecular Biosciences (GGNB). NB Paul is a member of the Molecular Medicine PhD program, part of the Göttingen Graduate Center for Neurosciences, Biophysics, and Molecular Biosciences (GGNB). T Beißbarth and M Haubrock are members of the Göttingen Campus Institute Data Science (CIDAS). Also, we thank the ENCODE Consortium and the ENCODE production laboratories for generating the particular datasets. The research was funded by Deutsche Forschungsgemeinschaft (DFG) (KFO5002 to M Sahrhage and SFB1002 to NB Paul). We acknowledge support by the Open Access Publication Funds/transformative agreements of the Göttingen University.

Author Contributions

M Sahrhage: conceptualization, resources, data curation, software, formal analysis, validation, investigation, visualization, methodology, project administration, and writing—original draft, review, and editing.

NB Paul: conceptualization, software, and writing—review and editing.

T Beißbarth: resources, supervision, funding acquisition, methodology, and project administration.

M Haubrock: conceptualization, data curation, software, supervision, methodology, project administration, and writing—review and editing.

Conflict of Interest Statement

The authors declare that they have no conflict of interest.

* Received September 18, 2023.

* Revision received May 15, 2024.

* Accepted May 16, 2024.

* © 2024 Sahrhage et al.

References

Kornberg RD (1974) Chromatin structure: A repeating unit of histones and DNA. Science 184: 868–871. doi:10.1126/science.184.4139.868

Luger K, Mäder AW, Richmond RK, Sargent DF, Richmond TJ (1997) Crystal structure of the nucleosome core particle at 2.8 Å resolution. Nature 389: 251–260. doi:10.1038/38444

Nora EP, Lajoie BR, Schulz EG, Giorgetti L, Okamoto I, Servant N, Piolot T, van Berkum NL, Meisig J, Sedat J, et al. (2012) Spatial partitioning of the regulatory landscape of the X-inactivation centre. Nature 485: 381–385. doi:10.1038/nature11049

Dixon JR, Selvaraj S, Yue F, Kim A, Li Y, Shen Y, Hu M, Liu JS, Ren B (2012) Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature 485: 376–380. doi:10.1038/nature11082

Clarkson CT, Deeks EA, Samarista R, Mamayusupova H, Zhurkin VB, Teif VB (2019) CTCF-dependent chromatin boundaries formed by asymmetric nucleosome arrays with decreased linker length. Nucleic Acids Res 47: 11181–11196. doi:10.1093/nar/gkz908

Zaret KS, Carroll JS (2011) Pioneer transcription factors: Establishing competence for gene expression. Genes Dev 25: 2227–2241. doi:10.1101/gad.176826.111

Slattery M, Zhou T, Yang L, Dantas Machado AC, Gordân R, Rohs R (2014) Absence of a simple code: How transcription factors read the genome. Trends Biochem Sci 39: 381–399. doi:10.1016/j.tibs.2014.07.002

Jiang C, Pugh BF (2009) Nucleosome positioning and gene regulation: Advances through genomics. Nat Rev Genet 10: 161–172. doi:10.1038/nrg2522

Lai WKM, Pugh BF (2017) Understanding nucleosome dynamics and their links to gene expression and DNA replication. Nat Rev Mol Cell Biol 18: 548–562. doi:10.1038/nrm.2017.47

Oberbeckmann E, Niebauer V, Watanabe S, Farnung L, Moldt M, Schmid A, Cramer P, Peterson CL, Eustermann S, Hopfner KP, et al. (2021) Ruler elements in chromatin remodelers set nucleosome array spacing and phasing. Nat Commun 12: 3232. doi:10.1038/s41467-021-23015-0

Kornblihtt AR, Schor IE, Alló M, Dujardin G, Petrillo E, Muñoz MJ (2013) Alternative splicing: A pivotal step between eukaryotic transcription and translation. Nat Rev Mol Cell Biol 14: 153–165. doi:10.1038/nrm3525

Fong N, Kim H, Zhou Y, Ji X, Qiu J, Saldi T, Diener K, Jones K, Fu XD, Bentley DL (2014) Pre-mRNA splicing is facilitated by an optimal RNA polymerase II elongation rate. Genes Dev 28: 2663–2676. doi:10.1101/gad.252106.114

Huang H, Chen J, Liu H, Sun X (2013) The nucleosome regulates the usage of polyadenylation sites in the human genome. BMC Genomics 14: 912. doi:10.1186/1471-2164-14-912

Barski A, Cuddapah S, Cui K, Roh TY, Schones DE, Wang Z, Wei G, Chepelev I, Zhao K (2007) High-resolution profiling of histone methylations in the human genome. Cell 129: 823–837. doi:10.1016/j.cell.2007.05.009

Bannister AJ, Kouzarides T (2011) Regulation of chromatin by histone modifications. Cell Res 21: 381–395. doi:10.1038/cr.2011.22

Struhl K, Segal E (2013) Determinants of nucleosome positioning. Nat Struct Mol Biol 20: 267–273. doi:10.1038/nsmb.2506

Zhang Y, Moqtaderi Z, Rattner BP, Euskirchen G, Snyder M, Kadonaga JT, Liu XS, Struhl K (2009) Intrinsic histone-DNA interactions are not the major determinant of nucleosome positions in vivo. Nat Struct Mol Biol 16: 847–852. doi:10.1038/nsmb.1636

Kaplan N, Moore IK, Fondufe-Mittendorf Y, Gossett AJ, Tillo D, Field Y, LeProust EM, Hughes TR, Lieb JD, Widom J, et al. (2009) The DNA-encoded nucleosome organization of a eukaryotic genome. Nature 458: 362–366. doi:10.1038/nature07667

Satchwell SC, Drew HR, Travers AA (1986) Sequence periodicities in chicken nucleosome core DNA. J Mol Biol 191: 659–675. doi:10.1016/0022-2836(86)90452-3

Segal E, Fondufe-Mittendorf Y, Chen L, Thåström A, Field Y, Moore IK, Wang JPZ, Widom J (2006) A genomic code for nucleosome positioning. Nature 442: 772–778. doi:10.1038/nature04979

Gaffney DJ, McVicker G, Pai AA, Fondufe-Mittendorf YN, Lewellen N, Michelini K, Widom J, Gilad Y, Pritchard JK (2012) Controls of nucleosome positioning in the human genome. PLoS Genet 8: e1003036. doi:10.1371/journal.pgen.1003036

Rohs R, West SM, Sosinsky A, Liu P, Mann RS, Honig B (2009) The role of DNA shape in protein–DNA recognition. Nature 461: 1248–1253. doi:10.1038/nature08473

Zhou T, Yang L, Lu Y, Dror I, Dantas Machado AC, Ghane T, Di Felice R, Rohs R (2013) DNAshape: A method for the high-throughput prediction of DNA structural features on a genomic scale. Nucleic Acids Res 41: W56–W62. doi:10.1093/nar/gkt437

Morozov AV, Fortney K, Gaykalova DA, Studitsky VM, Widom J, Siggia ED (2009) Using DNA mechanics to predict in vitro nucleosome positions and formation energies. Nucleic Acids Res 37: 4707–4722. doi:10.1093/nar/gkp475

Basu A, Bobrovnikov DG, Cieza B, Arcon JP, Qureshi Z, Orozco M, Ha T (2022) Deciphering the mechanical code of the genome and epigenome. Nat Struct Mol Biol 29: 1178–1187. doi:10.1038/s41594-022-00877-6

Segal E, Widom J (2009) Poly(dA:dT) tracts: Major determinants of nucleosome organization. Curr Opin Struct Biol 19: 65–71. doi:10.1016/j.sbi.2009.01.004

Liu H, Zhang R, Xiong W, Guan J, Zhuang Z, Zhou S (2014) A comparative evaluation on prediction methods of nucleosome positioning. Brief Bioinform 15: 1014–1027. doi:10.1093/bib/bbt062

Teif VB (2016) Nucleosome positioning: Resources and tools online. Brief Bioinform 17: 745–757. doi:10.1093/bib/bbv086

Han G-S, Li Q, Li Y (2021) Comparative analysis and prediction of nucleosome positioning using integrative feature representation and machine learning algorithms. BMC Bioinformatics 22: 129. doi:10.1186/s12859-021-04006-w

Shtumpf M, Piroeva KV, Agrawal SP, Jacob DR, Teif VB (2022) NucPosDB: A database of nucleosome positioning in vivo and nucleosomics of cell-free DNA. Chromosoma 131: 19–28. doi:10.1007/s00412-021-00766-9

Xi L, Fondufe-Mittendorf Y, Xia L, Flatow J, Widom J, Wang J-P (2010) Predicting nucleosome positioning using a duration Hidden Markov Model. BMC Bioinformatics 11: 346. doi:10.1186/1471-2105-11-346

Fang K, Li T, Huang Y, Jin VX (2021) NucHMM: A method for quantitative modeling of nucleosome organization identifying functional nucleosome states distinctly associated with splicing potentiality. Genome Biol 22: 250. doi:10.1186/s13059-021-02465-1

Guo S-H, Deng EZ, Xu LQ, Ding H, Lin H, Chen W, Chou KC (2014) iNuc-PseKNC: A sequence-based predictor for predicting nucleosome positioning in genomes with pseudo k-tuple nucleotide composition. Bioinformatics 30: 1522–1529. doi:10.1093/bioinformatics/btu083

Di Gangi M, Lo Bosco G, Rizzo R (2018) Deep learning architectures for prediction of nucleosome positioning from sequences data. BMC Bioinformatics 19: 418. doi:10.1186/s12859-018-2386-9

Zhang J, Peng W, Wang L (2018) LeNup: Learning nucleosome positioning from DNA sequences with improved convolutional neural networks. Bioinformatics 34: 1705–1712. doi:10.1093/bioinformatics/bty003

Routhier E, Pierre E, Khodabandelou G, Mozziconacci J (2021) Genome-wide prediction of DNA mutation effect on nucleosome positions for yeast synthetic genomics. Genome Res 31: 317–326. doi:10.1101/gr.264416.120

Han G-S, Li Q, Li Y (2022) Nucleosome positioning based on DNA sequence embedding and deep learning. BMC Genomics 23: 301. doi:10.1186/s12864-022-08508-6

Tsompana M, Buck MJ (2014) Chromatin accessibility: A window into the genome. Epigenetics Chromatin 7: 33. doi:10.1186/1756-8935-7-33

Schones DE, Cui K, Cuddapah S, Roh TY, Barski A, Wang Z, Wei G, Zhao K (2008) Dynamic regulation of nucleosome positioning in the human genome. Cell 132: 887–898. doi:10.1016/j.cell.2008.02.022

Boyle AP, Davis S, Shulha HP, Meltzer P, Margulies EH, Weng Z, Furey TS, Crawford GE (2008) High-resolution mapping and characterization of open chromatin across the genome. Cell 132: 311–322. doi:10.1016/j.cell.2007.12.014

ENCODE Project Consortium (2012) An integrated encyclopedia of DNA elements in the human genome. Nature 489: 57–74. doi:10.1038/nature11247

Breiman L (2001) Random forests. Machine Learn 45: 5–32. doi:10.1023/a:1010933404324

Chiu T-P, Comoglio F, Zhou T, Yang L, Paro R, Rohs R (2016) DNAshapeR: An R/bioconductor package for DNA shape prediction and feature encoding. Bioinformatics 32: 1211–1213. doi:10.1093/bioinformatics/btv735

Ernst J, Kellis M (2012) ChromHMM: Automating chromatin-state discovery and characterization. Nat Methods 9: 215–216. doi:10.1038/nmeth.1906

Karolchik D, Hinrichs AS, Furey TS, Roskin KM, Sugnet CW, Haussler D, Kent WJ (2004) The UCSC Table Browser data retrieval tool. Nucleic Acids Res 32: 493D–D496. doi:10.1093/nar/gkh103

Choi K-J, Quan MD, Qi C, Lee JH, Tsoi PS, Zahabiyon M, Bajic A, Hu L, Prasad BVV, Liao SCJ, et al. (2022) NANOG prion-like assembly mediates DNA bridging to facilitate chromatin reorganization and activation of pluripotency. Nat Cell Biol 24: 737–747. doi:10.1038/s41556-022-00896-x

Aubrey BJ, Strasser A, Kelly GL (2016) Tumor-suppressor functions of the TP53 pathway. Cold Spring Harb Perspect Med 6: a026062. doi:10.1101/cshperspect.a026062

Tanaka H, Takizawa Y, Takaku M, Kato D, Kumagawa Y, Grimm SA, Wade PA, Kurumizaka H (2020) Interaction of the pioneer transcription factor GATA3 with nucleosomes. Nat Commun 11: 4136. doi:10.1038/s41467-020-17959-y

Takaku M, Grimm SA, Shimbo T, Perera L, Menafra R, Stunnenberg HG, Archer TK, Machida S, Kurumizaka H, Wade PA (2016) GATA3-dependent cellular reprogramming requires activation-domain dependent recruitment of a chromatin remodeler. Genome Biol 17: 36. doi:10.1186/s13059-016-0897-0

Owens N, Papadopoulou T, Festuccia N, Tachtsidi A, Gonzalez I, Dubois A, Vandormael-Pournin S, Nora EP, Bruneau BG, Cohen-Tannoudji M, et al. (2019) CTCF confers local nucleosome resiliency after DNA replication and during mitosis. Elife 8: e47898. doi:10.7554/eLife.47898

Bai L, Morozov AV (2010) Gene regulation by nucleosome positioning. Trends Genet 26: 476–483. doi:10.1016/j.tig.2010.08.003

Herrmann CJ, Schmidt R, Kanitz A, Artimo P, Gruber AJ, Zavolan M (2020) PolyASite 2.0: A consolidated atlas of polyadenylation sites from 3′ end sequencing. Nucleic Acids Res 48: D174–D179. doi:10.1093/nar/gkz918

Samee MAH, Bruneau BG, Pollard KS (2019) A de novo shape motif discovery algorithm reveals preferences of transcription factors for DNA shape beyond sequence motifs. Cell Syst 8: 27–42.e6. doi:10.1016/j.cels.2018.12.001

Yu G, Wang L, He QY (2015) ChIPseeker: An R/bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics 31: 2382–2383. doi:10.1093/bioinformatics/btv145

Tompitak M, Vaillant C, Schiessel H (2017) Genomes of multicellular organisms have evolved to attract nucleosomes to promoter regions. Biophys J 112: 505–511. doi:10.1016/j.bpj.2016.12.041

Klemm SL, Shipony Z, Greenleaf WJ (2019) Chromatin accessibility and the regulatory epigenome. Nat Rev Genet 20: 207–220. doi:10.1038/s41576-018-0089-8

Tillo D, Kaplan N, Moore IK, Fondufe-Mittendorf Y, Gossett AJ, Field Y, Lieb JD, Widom J, Segal E, Hughes TR (2010) High nucleosome occupancy is encoded at human regulatory sequences. PLoS One 5: e9129. doi:10.1371/journal.pone.0009129

Barozzi I, Simonatto M, Bonifacio S, Yang L, Rohs R, Ghisletti S, Natoli G (2014) Coregulation of transcription factor binding and nucleosome occupancy through DNA features of mammalian enhancers. Mol Cell 54: 844–857. doi:10.1016/j.molcel.2014.04.006

Zhu F, Farnung L, Kaasinen E, Sahu B, Yin Y, Wei B, Dodonova SO, Nitta KR, Morgunova E, Taipale M, et al. (2018) The interaction landscape between transcription factors and the nucleosome. Nature 562: 76–81. doi:10.1038/s41586-018-0549-5

He HH, Meyer CA, Shin H, Bailey ST, Wei G, Wang Q, Zhang Y, Xu K, Ni M, Lupien M, et al. (2010) Nucleosome dynamics define transcriptional enhancers. Nat Genet 42: 343–347. doi:10.1038/ng.545

Andersson R, Enroth S, Rada-Iglesias A, Wadelius C, Komorowski J (2009) Nucleosomes are well positioned in exons and carry characteristic histone modifications. Genome Res 19: 1732–1741. doi:10.1101/gr.092353.109

Kogan S, Trifonov EN (2005) Gene splice sites correlate with nucleosome positions. Gene 352: 57–62. doi:10.1016/j.gene.2005.03.004

Chodavarapu RK, Feng S, Bernatavichute YV, Chen PY, Stroud H, Yu Y, Hetzel JA, Kuo F, Kim J, Cokus SJ, et al. (2010) Relationship between nucleosome positioning and DNA methylation. Nature 466: 388–392. doi:10.1038/nature09147

Mavrich TN, Ioshikhes IP, Venters BJ, Jiang C, Tomsho LP, Qi J, Schuster SC, Albert I, Pugh BF (2008) A barrier nucleosome model for statistical positioning of nucleosomes throughout the yeast genome. Genome Res 18: 1073–1083. doi:10.1101/gr.078261.108

Snyder M, Raha D, Johnson S, Winters E, Sidow A, Weng Z, Smith C, Lacroute P, Cayting P, Kundaje A (2011) General Protocol - track description for UCSC Genome Browser composite track hg19/wgEncodeSydhNsome. Available at: https://www.encodeproject.org/documents/912dc886-d576-48ac-aa75-58ddaee6d9ff/@@download/attachment/wgEncodeSydhNsome.html.pdf [accessed: 01 September 2023].

Boyle A, Song L, Crawford G, Giresi P, Lieb J, Furey D, Liu Z, McDaniell R, Lee B-K, Iyer V, et al. (2011) General Protocol - track description for UCSC Genome Browser composite track hg19/wgEncodeOpenChromDnase. Available at: https://www.encodeproject.org/documents/6878f54f-b81e-4064-a6ae-f78c76ed7f85/@@download/attachment/wgEncodeOpenChromDnase.release3.html.pdf [accessed: 01 September 2023].

Onate K (2016) Pipeline Protocol - ChIP pipeline: Details for both histone and transcription factor. Available at: https://www.encodeproject.org/documents/7009beb8-340b-4e71-b9db-53bb020c7fe2/@@download/attachment/ChIP-seq_pipeline_overview.pdf [accessed: 01 September 2023].

Pollard KS, Hubisz MJ, Rosenbloom KR, Siepel A (2010) Detection of nonneutral substitution rates on mammalian phylogenies. Genome Res 20: 110–121. doi:10.1101/gr.097857.109

O’Leary NA, Wright MW, Brister JR, Ciufo S, Haddad D, McVeigh R, Rajput B, Robbertse B, Smith-White B, Ako-Adjei D, et al. (2016) Reference sequence (RefSeq) database at NCBI: Current status, taxonomic expansion, and functional annotation. Nucleic Acids Res 44: D733–D745. doi:10.1093/nar/gkv1189

Lawrence M, Gentleman R, Carey V (2009) rtracklayer: An R package for interfacing with genome browsers. Bioinformatics 25: 1841–1842. doi:10.1093/bioinformatics/btp328

R Core Team (2023) R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing.

Wickham H (2016) Ggplot2: Elegant graphics for data analysis. In Use R!, 2nd ed. Cham, Switzerland: Springer International Publishing.

Lawrence M, Huber W, Pagès H, Aboyoun P, Carlson M, Gentleman R, Morgan MT, Carey VJ (2013) Software for computing and annotating genomic ranges. PLoS Comput Biol 9: e1003118. doi:10.1371/journal.pcbi.1003118

Sing T, Sander O, Beerenwinkel N, Lengauer T (2005) ROCR: Visualizing classifier performance in R. Bioinformatics 21: 3940–3941. doi:10.1093/bioinformatics/bti623

Liaw A, Wiener M (2002) Classification and regression by randomForest. R News 2: 18–22.

Amato D, Bosco GL, Rizzo R (2020) CORENup: A combination of convolutional and recurrent deep neural networks for nucleosome positioning identification. BMC Bioinformatics 21: 326. doi:10.1186/s12859-020-03627-x

Valouev A, Johnson SM, Boyd SD, Smith CL, Fire AZ, Sidow A (2011) Determinants of nucleosome organization in primary human cells. Nature 474: 516–520. doi:10.1038/nature10002

Liao Y, Smyth GK, Shi W (2019) The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. Nucleic Acids Res 47: e47. doi:10.1093/nar/gkz114