Free energy-based model of CTCF-mediated chromatin looping in the human genome

Highlights

• Dynamic programming algorithm systematically searches chromatin free energy landscape.

• Free energy minimization generates folded optimal and suboptimal 2D meta-structures.

• 3D structures are generated from 2D meta-structure including multipair interactions.

• The method identifies the contacts that are significant in the 2D contact map.

• CTCF contacts and free energy correlate strongly with regions containing euchromatin.

Abstract

In recent years, high-throughput techniques have revealed considerable structural organization of the human genome with diverse regions of the chromatin interacting with each other in the form of loops. Some of these loops are quite complex and may encompass regions comprised of many interacting chain segments around a central locus. Popular techniques for extracting this information are chromatin interaction analysis by paired-end tag sequencing (ChIA-PET) and high-throughput chromosome conformation capture (Hi-C). Here, we introduce a physics-based method to predict the three-dimensional structure of chromatin from population-averaged ChIA-PET data. The approach uses experimentally-validated data from human B-lymphoblastoid cells to generate 2D meta-structures of chromatin using a dynamic programming algorithm that explores the chromatin free energy landscape. By generating both optimal and suboptimal meta-structures we can calculate both the free energy and additionally the relative thermodynamic probability. A 3D structure prediction program with applied restraints then can be used to generate the tertiary structures. The main advantage of this approach for population-averaged experimental data is that it provides a way to distinguish between the principal and the spurious contacts. This study also finds that euchromatin appear to have rather precisely regulated 2D meta-structures compared to heterochromatin. The program source-code is available at https://github.com/plewczynski/looper.

Introduction

Recent high-throughput sequencing methods on chromatin fibers have revealed considerable hierarchical organization and packaging of the chromatin in the form of loop-like structures in humans [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], other vertebrates such as mice [14], or zebrafish [15], insects like Drosophila melanogaster [16], and even plants [17]. The 3D topology of chromatin not only helps to establish the volume of the nucleus, it is thought to have a functional role; e.g. in regulating gene expression by modulating the spatial distance between promoters and enhancers [2], [3], [4], [18], [19], [20], [21]. Alterations in the arrangement of some of the loops have been linked to diseases like cancer [22] and autoimmune diseases [23].

A key feature in humans and other vertebrates is that the DNA loops arise as a consequence of protein-mediated binding of distinct DNA regions. A pivotal role in this process is played by CTCF: an 11 zinc-fingers protein that binds to the CCCTC motif [2], [18], [21], [24]. Proteins, which are often found concomitant with the CTCF protein around the DNA loops, include cohesin, Yin Yang 1, the RNA polymerase II complex (RNAPII), or CCCTC-binding factor like protein (CTCFL: also called BORIS) [2], [18], [21], [24], [25]. The latter is expressed in germ cells and some cancers and mediates the formation of a subset of cell-specific loops [26]. YY1-mediated loops bring promoters and enhancers into close spatial proximity [27]. RNAPII also appears to be strongly associated with the presence and formation of chromatin loops; though often shorter and more tissue specific [2]. How the loops are formed is an open question. In some works, it was shown that CTCF forms a dimer that can bring together two distal regions of DNA [18]. Another possible mechanism is the loop extrusion model, which proposes that the pair of CTCFs serves as a topological obstacle for the cohesin ring traversing along the DNA fiber [28], [29]. RNAPII may also play an important role in the active translocation of the cohesin.

There are multiple methods that have been developed to obtain genome-wide maps of chromatin organization [30], [31]. The most common approach is high-throughput chromosome conformation capture (Hi-C), which identifies all types of interactions between diverse parts of the chromatin chain. Based on Hi-C data, the genome can be divided into substructures known as topologically associated domains – TADs [3], [6], [13], [32]. A TAD is essentially contiguous regions of the chromatin in which loci tend to interact far more frequently with each other than with loci outside the region. Another method is chromatin interaction analysis by paired-end tag sequencing (ChIA-PET), which is designed to measure only the interactions mediated by specific protein factors like CTCF or RNAPII [33], [34]. This is achieved by introducing the chromatin immunoprecipitation (ChIP) step that targets specific proteins associated with the protein-DNA complexes [33]. Furthermore, ChIA-PET permits identifying specific chromatin loops at a resolution of DNA base pairs, in principle. The strength of these loops is expressed in terms of the number of PET counts; i.e., the number of unique pairs of reads confirming that two regions are joined together. Usually only loops above a specific PET count threshold are reported; however, even loops with PET count equal to 1 (called ‘singletons’), which may arise due to transient interactions with the chromatin chain, are sometimes used for example to propose chromatin 3D structure [1]. Two interacting loci that form a given loop are called anchors [2], [31]. The length (or “span”) of the ChIA-PET-proposed loops can vary significantly; between 8 kilobase-pairs (kbp) up to 1 Mbp. Similarly to Hi-C, with ChIA-PET, the genome can be partitioned into segments called chromatin contact domains (CCDs). A CCD is a genome fragment that is continously covered by interconneted loops and, for the GM12878 cell line, 2,267 CCDs were identified. The ChIA-PET experiments on the GM12878 cell line produced 31 million reads [2]. The data is usually binned into segments of 1- to 5-kbp. In general, a resolution of 5 kbp produces a reasonable 2D contact map (or heat map). Beside the two abovementioned methods, several other techniques have been introduced recently to analyze 3D genome organization like ChIA-Drop [35], GAM [36], or PLAC-Seq [37]. Some of these strive to report the chromatin structure from individual cells rather than a population-averaged picture observed with older techniques.

Chromatin loops can be described in terms of motifs. For example, Fig. 1A depicts a pair of CTCFs interacting with cohesin, a ring-shaped multidomain structural protein, resulting in a simple loop. As the CTCF binding motifs in both anchors are in a parallel orientation, we call this loop “convergent”. This type of loop is the most frequently observed substructure in the human genome. Based on the PET count, it also can be considered the strongest of interactions. “Divergent” loops occur when the CTCF binding motifs adopt an antiparallel orientation, thus influencing the relative orientation of both CTCFs bound to this DNA fragment (Fig. 1B). When the direction of the motifs is parallel, two types of structures of equal tendency and intermediate strength are suggested: tandem right (Fig. 1C) and tandem left (Fig. 1D). The CTCF ChIA-PET data for the GM12878 cell line indicate that around 65% of the loops are convergent, roughly 33% are tandem and remaining divergent loops amount to only around 2%. Similar values were obtained for Hi-C data for the same cell line [2], [21]. CTCF structures also appear to group into islands (also called multipair interactions, as shown in Fig. 1E) when more than two chromatin chains interact at a single spatial locus. In this case the central region is the largely insoluble and inactive part (heterochromatin) whereas the regions jutting out are more accessible to transcription factors and therefore active (euchromatin).

In recent years, there has been a considerable interest in finding ways to physically model the distribution of chromatin structure using various approaches such as population-based analysis [38], polymer based models using molecular dynamics (MD) simulations [2], [11], [12], [39], [40], or Monte Carlo (MC) simulation techniques [1], [8], and some thermodynamics based methods [16], [41], [42]. These models provide 3D structures of the chromatin based on the population-averaged experimental data. However, MD or MC simulation techniques do not guarantee an exhaustive search of such a large landscape as is observed for chromatin, nor is it assured that a conformation found in this way is the most stable structure [38]. More recently, a maximum entropy method has been introduced that explores the phase space of the polymer chain under the constraints of pair-wise interactions and uses these contacts to fit a proposed 3D potential using maximum entropy [43], [44].

Here, we introduce a different approach where we use thermodynamics on 2D meta-structures to produce ranked, optimal and suboptimal structures. In our model, we assume that the number of observed PET counts associated with a particular contact point is a measure of the likelihood of that interaction. Given so, it should be a function of the relative stability of the DNA loops and can be represented as a binding free energy. For finite, uniquely definable interactions, a thermodynamic potential can be developed to estimate the likelihood of a given 2D meta-structure based upon the free energy of each motif. These attractive interactions were combined with the repulsive effects caused by entropy loss due to folding and the formation of contacts; a concept originally introduced as the cross linking entropy ([45], [46], [47], [48] and references therein). This approach provides a list of optimal and suboptimal 2D meta-structures from the collection of possible loops. To systematically rank various combinations of loops, we employed a dynamic programming algorithm (DPA) strategy to optimize the arrangement of these pair contacts. Using heuristic adaptations to the DPA, the method is designed to handle complex motifs such as pseudoknots involving both parallel and antiparallel chain interactions (see Section 2.3). Our model works largely from the opposite direction of what is proposed by Di Pierro et al. [43], which is a top-down approach. Combining the statistical potential with the entropy loss, this model provides a quantitative, bottom-up integration procedure that is consistent with the Di Pierro et al. maximum entropy model.

Unlike many protein and RNA structures where there is a very well-defined native 3D structure; chromatin conformations tend to be far more plastic. These highly flexible loops effectively confine the enclosed chromatin chain to a subset of possible interactions; increasing the probability that some of the observed contacts are spurious interactions mixed with the genuine contacts. During the cross-linking process, commonly used in experimental techniques, the chromatin is fixed in a specific 3D conformation. Consequently, CTCF-cohesin complexes in transient locations on the chromatin chain will end up fused at those locations and pulled down with the actual, biologically relevant contacts. If spurious contacts are present in the population averaged data, then they will also appear in single-cell data and the priority of the contacts cannot be clearly established. This presents some challenges in interpreting 2D contact maps. Our method provides a way to “filter” and rank the population-averaged and single-cell data into a set of human readable optimal and suboptimal structures and generates 3D structures and trajectories from that information. This modeling approach tends to be consistent with the observed behavior of DNA. The 3D structures are not calculated from rigid harmonic potentials, rather, the structures are systematically ranked using the DPA and then the 3D structure is constructed using restraints associated with these predicted contacts. The thermodynamic approach also tends to guarantee that the loss of pairing probability with genomic distance n is consistent with the observed behavior: p(n) ∝ 1/n^α, where α = 1 would represent the fractal globule case with respect to the volume of the nucleus [6].

In this approach, we also focus on representing chromatin motifs as 2D meta-structures. This representation approach is helpful both in visualizing chromatin structures, and in comparing differences between two related genomes. This helps focus any 3D visualizations on the relevant changes in a given genome. Furthermore, computation with our approach generates 1D structural representations, which are machine readable, and can be automated to generate 3D conformations.

In summary we propose a novel method to help filter out spurious contacts obtain with methods designed to analyze genome-wide chromatin organization and to obtain the optimal structures of chromatin. As proof of principle, we focus on population-averaged ChIA-PET data from the GM12878 cell line [1], [2]. We applied the free energy prediction approach to obtain optimal and suboptimal meta-structures that reveal an ensemble of chromatin structural motifs. To validate our approach, we predicted both the optimal structure and all the other major suboptimal meta-structures of chromatin regions for the publicly available list of chromatin contact domains (CCDs), as well as for the members of A and B compartments. We observe that, for segments of chromatin that are active, there are more suboptimal structures that tend to have a more diverse distribution of free energies than regions that are inactive. Thus, with this approach we can correctly recognize regions of chromatin with high mobility identified in independent experimental methods.