Mapping the spectrum of 3D communities in human chromosome conformation capture data.

Several experiments show that the three dimensional (3D) organization of chromosomes affects genetic processes such as transcription and gene regulation. To better understand this connection, researchers developed the Hi-C method that is able to detect the pairwise physical contacts of all chromosomal loci. The Hi-C data show that chromosomes are composed of 3D compartments that range over a variety of scales. However, it is challenging to systematically detect these cross-scale structures. Most studies have therefore designed methods for specific scales to study foremost topologically associated domains (TADs) and A/B compartments. To go beyond this limitation, we tailor a network community detection method that finds communities in compact fractal globule polymer systems. Our method allows us to continuously scan through all scales with a single resolution parameter. We found: (i) polymer segments belonging to the same 3D community do not have to be in consecutive order along the polymer chain. In other words, several TADs may belong to the same 3D community. (ii) CTCF proteins-a loop-stabilizing protein that is ascribed a big role in TAD formation-are well correlated with community borders only at one level of organization. (iii) TADs and A/B compartments are traditionally treated as two weakly related 3D structures and detected with different algorithms. With our method, we detect both by simply adjusting the resolution parameter. We therefore argue that they represent two specific levels of a continuous spectrum 3D communities, rather than seeing them as different structural entities.

Introduction

Experiments that detect the physical contacts between DNA loci in the nucleus, such as Hi-C[1,2], show that DNA in the nucleus is not a randomly folded polymer. Rather, across cell types and organisms, Hi-C experiments reveal that chromosomes are built up by a network of three dimensional (3D) compartments.

At the mega (106)-base-pair scale, two types of coexisting structures stand out. In one, all chromosome loci belong to one of the two so-called A/B compartments, where the chromatin in one compartment is generally more open, accessible, and actively transcribed than the other. In the second type, linear subsections of the genome assemble into topological domains[3,4], often referred to as topologically associated domains (TADs)[2-4]. Plotting Hi-C data as a heat map, TADs show up as local regions with sharp borders with more internal than external contacts. The positions of these borders are correlated with several genetic processes, such as transcription, localization of some epigenetic marks, and DNA-binding profiles of several proteins–most notably CTCF and cohesin[3,4].

The methods that detect TADs are not the same as those that find A/B compartments. Therefore, TADs and A/B compartments are treated as different 3D structures that are only weakly related to each other. Just as for TADs, there are several algorithms tailored for detecting A/B compartments[1,5], each with their strengths and weaknesses.

To algorithmically detect TAD-like structures, there exists by now a menagerie of network[6-10] and clustering approaches[11-17]. Arguably these methods yield overlapping results, but it is unclear by how much. In particular, some methods cannot deal with TAD-within-TAD hierarchies that become apparent when zooming in TADs in highly resolved Hi-C maps. This means that there is not a universal definition for what a TAD really is.

Some network approaches are based on community detection methods that are related to what we use here. In Cabreros et al.[7], the authors suggest a method based on the stochastic block model[16,17], which is another side of the network community detection field compared to the modularity maximization method that we use here (but as shown recently[18], they are connected). However, there are some limitations in their approach. For example, they binarize the Hi-C data (‘no connection’ or ‘connection’ with a rather arbitrarily chosen thresholds) thereby discarding contact frequency variations in the Hi-C data. Furthermore, there is no comparative study connecting their communities to biological factors or any mechanistic models. Wang et al.[8] takes a step in this direction, but the method to detect the TADs itself relies on biological factor data and nontrivial threshold criteria.

To overcome some of these problems, we start by acknowledging that the chromosomes have a richer 3D organization than simply TADs and A/B compartments. These are just two examples. To capture this, we develop network-based method that allows us to scan through 3D structures on all scales with a resolution parameter. In particular, our approach is based on the GenLouvain method, originally designed for network community detection. For a specific value of the resolution parameter, the method finds the optimal community structure with respect to a null model of the network that has to be specified beforehand. Based on the physics of compact polymer globules, we put forward a null model that is consistent with the average contact probabilities in real Hi-C data[1]. This goes beyond previous Louvain-like studies[14,15] that treat the Hi-C data as a network with random connections.

Furthermore, most studies, such as Yan et al. and Norton et al.[14,15], treat TADs as linear contiguous sequences of chromatin. This restriction overrides the GenLouvain algorithm’s ability to find the (not necessarily contiguous) optimal community structure in the data set[19,20]. We remove this restriction in our study. Therefore, to reduce confusion, we will not use the term TAD, but rather the 3D community for the cluster of nodes that comes out of the GenLouvain algorithm since they are not necessarily linear contiguous sequences.

Methods

We use the GenLouvain algorithm[21] to detect communities in the Hi-C maps (this is an extension of the original Louvain method[22]). This algorithm offers the possibility to find communities at several scales with a single resolution parameter γ. Contrasting other methods with similar features, the spectrum of communities that we detect is not necessarily hierarchical or nested, as in e.g. Wang et al., Fraser et al., and Bianco et al.[8,23,24]. Instead, two different values of γ give two different collections of communities, and these do not necessarily have anything to do with each other.

To find the communities, the GenLouvain method tries putting the nodes into different communities to maximize the so-called modularity function Q. This function quantifies by how much more dense the connections are within the communities compared to what they would be in a particular type of network, say a random network. GenLouvain’s modularity function iswhere the sum runs over all nodes in the network, in our case genomic loci, and A is the number of contacts between all node pairs i and j (A is the adjacency matrix in standard network terminology). The resolution parameter γ controls the overall scale of communities we would like to detect, where larger values of γ correspond to smaller communities. Furthermore, is the null-model term that is network-type-specific, and the difference therefore measures how strongly nodes i and j are connected in the real network, compared to how strongly we expect them to be given . Finally, δ(g, g) is the Kronecker delta that is unity only if nodes i and j belong to the same community (otherwise it is zero), and m is a normalization factor so that Q goes between −1 and +1.

One of the most popular choices for is the Newman-Girvan (NG) null-model term for a random network, [25]. For a unweighted network where A is either zero or one, k and k are the number of links for nodes i and j (). Simply put, the NG null-model assumes that the probability that i and j are connected is proportional to the product of their number of links. The same interpretation holds for weighted networks where k becomes the sum of weights on the edges connected to i (“strength” in network terminology).

However, the NG null-model is too rough an approximation to find communities in Hi-C maps, because it does not obey the well-established contact patterns that we know exist in DNA, or in fact any polymer system. For example, DNA is a long polymer where nodes are arranged in a linear sequence and then folded in 3D. This sets limitations for how frequently two pieces of DNA, or nodes, can join in 3D space (this is usually not a restriction in most networks). As was first discovered in Lieberman-Aiden et al.[1], and many papers thereafter[26,27], the contact probability between two nodes i and j decays as a power-law with the linear distance between them, that is . Based on this, we propose the following null-model:

The value of the decay exponent α is debatable, but we use α = 1. There are two main reasons. First, at the mega-base-pair scale of human DNA, the Hi-C data suggest that it is close to one[1] (α is also close to one in mice[27]). Second, in the next section we will study community detection in the fractal globule polymer (hence the superscript FG in ) where α = 1[28,29]. Nonetheless, we point out that our method does note rely on this choice, and α is in principle a free parameter.

3D communities in fractal globules

Before we investigate the human Hi-C data, we wish to better understand the types of 3D structures that our community detection method picks up. To do this, we use computer-generated fractal globule polymers (denoted by the “crumpled globule” in the original article[30]). This is a compact polymer that mimics the large scale structure of human chromosomes, in particular the scaling relations of the end-to-end distance and contact frequency[1]. The advantage of this approach is that we have the explicit 3D coordinates for every part of the polymer (because we made it), which allows us to visualize and analyze the 3D structure of the communities that we detect. For the Hi-C map, we only have pairwise contact frequencies.

Generating fractal globules with the conformation-dependent polymerization algorithm

As introduced in an article[31], there are several variants of fractal globules-like structures[2,32-34]. Due to simplicity and speed, we use the conformation-dependent polymerization (CDP) model[28]. In a nutshell, this is a Monte Carlo method that produces a fractal globule by simulating a biased random walk on lattice where the propagation probability depends on the entire walk’s trajectory over the lattice[28]. This yields on-lattice space-filling polymers. To generate off-lattice fractal globules, the structure is randomized with simulated annealing where the position of a monomer is randomly displaced under the constraint of a fixed inter-monomer distance. With properly chosen parameters, the CDP method produces a fractal globule with contact frequency (probability) that decays as ~s−1 (as it should[29]), where s is the contour distance along the polymer. To use the notation from Tamm et al.[28], we use the relation p ∝ (1 + An) with A = 104 for an unoccupied site and p to an occupied site is given by ε = 10−4. This result is consistent with the Hi-C data for the human genome at the mega-base-pair scale[1]. We show a realization of a CDP in Fig. 1(a).

Figure 1

3D communities in simulated fractal globules. (a, left) The 3D representation of fractal globule from the CDP algorithm. The colors highlight communities that we detect using Eq. (2), γ = 0.6. (a, right) The communities are not contiguous: the small globule is a subsection of the lower left part of polymer, and the stretched version shows the alternating communities (ABCAB). (b–d) The contact maps of simulated fractal globules after KR-normalization, with various resolution parameters: (b) γ = 0.4, (c) γ = 0.6, and (d) γ = 0.8. To show non-contiguous communities, we superimpose them as the squares; the same color indicates that they belong to the same 3D community. (e) The end-to-end distance for the fractal (FG, blue) and the equilibrium globules (EG, red) averaged over 200 polymer realizations. The triangles denote the end-to-end distance for community boundaries, and dashed lines represent the chain as a whole. To find the communities, we use γ = 0.4. The data are obtained from the simulation of 200 sample globules for each polymer model. The error bars show the standard error of the mean, and the two guided slopes (1/2 and 1/3) show the known scaling of equilibrium and fractal globules at intermediate length scales. (f) The same data as in the panel (e), where we scale the vertical axis with the radius of gyration where N is the total number of polymer segments and ri is the coordinate of the ith segment.

From every simulated CDP polymer, we construct a contact map by counting the number of contact events between polymer beads i and j. The contact refers to the case where the Euclidean distance between two beads are shorter than three lattice spacings. To obtain good statistics, we first generate an on-lattice polymer, and then during the annealing stage we register all contacts in 103 random variations of the on-lattice structure. In addition, we have confirmed that different threshold values for what we consider as a contact do not qualitatively alter our results presented below. Just as in the Hi-C experiments, in the final step we normalize the contact map with the KR-norm[35] so that each row and column sums to unity.

As a reference case, we use the equilibrium globule. This is a self-avoiding polymer in a closed spherical volume. When we generate equilibrium globules, we contain it in a volume with the same diameter as the fractal globules’ radius of gyration.

Structure of 3D communities in fractal globules

Figure 1(b) shows the contact map for one simulated CDP, where high pixel intensity indicates many contacts. Just as in real Hi-C maps, the simulated map contains locally concentrated contact domains along the diagonal, that is, along the polymer chain. To find these domains algorithmically, we put the contact map into our modified GenLouvain method. By varying the resolution parameter γ, we detect communities on various scales. On top of the contact map in Fig. 1(b–d), we overlay examples of 3D communities for γ = 0.4 (b), γ = 0.6 (c), and γ = 0.8 (d), where the boxes represent community boundaries. Based on these, we ask what the 3D structure of these communities is, and if and how they are different from the polymer as a whole.

Contrasting current views on TADs, we find that 3D communities do not have to be a contiguous polymer segment. Rather, linearly distant parts of the polymer can fold in 3D to form a community. We show this in Fig. 1(a) (left), where we mark the polymer segments belonging to the same community with the same color (γ = 0.6). In Fig. 1(a) (right), we cut out a subsection with three communities and stretch it out. Labeling the communities as A, B, and C, they are clearly ordered in a non-contiguous sequence: they appear as A-B-C-A-B rather than A-B-C.

Furthermore, because of the above-average contact frequencies inside a community, we would like to quantify how its 3D structure differs from the fractal globule polymer as a whole. To do this, we examine the scaling relation of the end-to-end distance–the Euclidean distance between the two boundary monomers defining that (contiguous) community–with respect to the subchain length s.

In Fig. 1(e), we show this relation for the community subchains (the blue triangles) and for all subchains (the blue dashed lines). It shows that the end-to-end distance for the entire globule grows as ~s1/3, as is expected for a space-filling curve (deviations for large-s comes from finite-size effects and insufficient statistics). For the communities, we notice that the end-to-end distances are systematically smaller than for a randomly chosen subchain. Our simulations even suggest that the scaling exponent is smaller than 1/3. Consistent properties are crosschecked in Fig. 1(f) where the end-to-end distance divided by the average chain size (2 × radius of gyration) is plotted. Overall, this shows that our method detects 3D communities that are compact substructures of the fractal globule. This observation supports that some TADs in chromosomes are end-closed loop structures[2]. For comparison, we made the same analysis for the equilibrium globule (EG). In Fig. 1(e), we see that the end-to-end distance for our 3D communities and all subchains have the same scaling: ~s1/2 for and ~s0 for , as we expect from an N monomer ideal chain. Concurrently, the end-to-end distance normalized by the radius of gyration (Fig. 1(f)) are almost same for both 3D communities and all subchains, showing a sharp contrast to the FG.

We further investigate the asphericity of the communities’ 3D structure in Supplementary Information (SI). We find that the community subchains in FGs tend to be more sphere-like than the FG itself, and the subchains in EGs. This leads to enhanced contact frequency between monomers within a community (Supplementary Fig. S1). Refer to SI for further discussion.

All these analyses support that our community detection method successfully identifies strongly interacting subchains in an FG from the contact map, and that these communities have polymeric properties that cannot be explained by the global expected behavior. Rather, they are more similar to features of TADs. Additionally, in the following section, we corroborate our method by comparing the communities in Hi-C maps detected by our method to the TAD data reported in the literature[2].

Community detection for the Hi-C map

Based on our community detection approach, we now proceed to analyze Hi-C data from human cells[2]. The data comes in the form of matrices where each entry represents the number contacts between two chromosome loci i and j. As is standard in the field, we normalized the data with the KR-norm[35] which balances the matrix such that every row and column sum to unity (we also used the KR-norm for the fractal globule contact data). The data is available in various resolutions, from 103 base pairs (1 kbp) to 106 base pairs (1 Mbp), but we used 100 kbp which is the scale where both TADs and A/B compartments can be detected[2].

In Fig. 2(a–c), we show that our algorithm detects differently sized contiguous blocks, or TADs, along the chromosome arms as we change the resolution parameter γ. Similar to the simulated data, it is clear that several TADs form 3D communities. To investigate how these TADs correspond to TADs defined in other studies, we compared the border locations of each contiguous block to TAD borders defined by Rao et al.[2] (the group that produced the dataset we use in this study) in Fig. 2(d). At γ ≈ 0.7, about 70% of the borders overlap. Moreover, several studies have shown that binding sites for the insulator protein CTCF are strongly correlated with TAD borders (e.g., in Dixon et al.[3]). We therefore check the overlap between CTCF binding sites (mapped by ENCODE[36]) and all borders at different γ values. We note that CTCF has the highest overlap also at γ ≈ 0.7 (Fig. 2(d)).

Figure 2

3D communities in real Hi-C data (chromosome 1). (a–c) Normalized Hi-C data with squares showing the structure of 3D communities. The black regions are the unmappable regions. The resolution parameter ranges from γ = 0.6 (a), γ = 0.7 (b), and γ = 0.8 (c). As in Fig. 1(b–d), we assign the same colors to those squares that belong to the same 3D community. It is clear that they are not contiguous sequences. (d) The fraction of community boundary points predicted by our method that coincide with the ones in Rao et al.[2] (the squares), and binding positions for CTCF (the circles) for different values of γ. (e) The average number of TADs for each community, and the number of communities as functions of γ. (f) The community division along chromosome 1 for three different values of γ. The purple squares represent the largest community in the panels (g)–(i), while the other colors indicate smaller communities. (g)–(i) Communities’ gene activity sorted by their relative size for different values of γ: (g) γ = 0.6, (h) γ = 0.65, and (i) γ = 0.7. The circles show the median RNA expression levels, and vertical lines are quartiles. We omit communities that are smaller than 50 nodes. We find the communities using Eqs (1 and 2).

As we discussed in the introduction, there are different algorithms that detect TADs. However, regardless of the definition that is used, the TADs that come out have substructures that we can interpret as TAD-within-TADs[12] with new sets of borders. Our algorithm lets researchers scan through all these TAD-within-TADs. According to Fig. 2(d), CTCF correlates well with TAD borders at γ ≈ 0.7. This level also coincides with TADs from Rao et al.[2]. Our method opens the possibility to investigate which family of proteins is important for different hierarchical levels.

We next investigate properties of the 3D communities. First, we count the number of communities and average number of TADs within each community at different γ (Fig. 2(e)) as well as their localization along the chromosome arms (Fig. 2(f)). The highest number of TADs/community is when γ ≈ 0.6 and γ ≈ 0.75 (about 17 TADs/community). Interestingly, for γ < 0.6, the chromosome consists of only two communities (but with several TADs within them).

Second, we map gene activity within each community using RNA-seq data from ENCODE[36] (GEO Acc. Nr: GSE88583). In Figs. 2(g–i) and Supplementary Figs. S2–S5, we show the average coverage of RNA-seq reads within communities for different γ. At small γ values where only two communities are defined, one community is clearly more active than the other. The active and less active communities are then split up into smaller communities as γ increases. Already in the original Hi-C paper by Lieberman-Aiden et al.[1], they used principal component analysis on the Hi-C data to partition the chromosomes into two classes. Denoting them by A/B compartments, they found that the A compartment contains transcriptionally active chromatin, whereas the B compartment is less active.

We took one step in this direction by looking at different types of chromatin. We use 15 chromatin states defined by ENCODE[37,38] to see what types of chromatin states are enriched in different communities at different γ values (Supplementary Fig. S6). These results show that at low values of γ one of the two large communities that are identified is enriched in chromatin states associated with transcription, while the other community is not. This again confirms that these represent the A and B compartments respectively. A very small third compartment is also visible that consists of centromere proximal repetitive regions. With increasing γ values, we can then see that first the active compartment (A) starts to split up into smaller compartments and then the less active parts of the genome (B) also start to split into smaller communities. It is clear that the genome does not split up into evenly sized sub-compartments, but rather the 3D space is dominated by a few large communities. It is striking and unprecedented that by tuning a single parameter, we detect both TADs and A/B compartments with the same algorithm. We therefore argue that TADs and A/B compartments are not two conceptually different organizational structures in the nucleus, but rather different ends of the same organizational spectrum.

Discussion

From Hi-C experiments, it is clear that inter-phase chromosomes are built up by a network of 3D compartments on various scales–from kilo(103)-base-pair sized loops to mega(106)-base-pair sized 3D structures. This pattern is consistent across organisms and cell types. To let researchers scan through the spectrum of 3D compartments, we have tailored the GenLouvain community detection method to find 3D communities in fractal globule polymer systems. Apart from verifying our method on computer-generated polymers, we have applied it to analyze human Hi-C data. First, we have found that chromatin segments belonging to the same 3D community do not have to be in next to each other along the DNA. In other words, several TADs can belong to the same 3D community. Second, we have found that CTCF proteins–a loop-stabilizing protein that is ascribed a big role in TAD formation–are only correlated well with community borders at one level of organization. It remains to find what other factors are important at higher or lower levels. Third, just by adjusting a single parameter (γ), our method picks up the two most prominent 3D compartments, TADs and A/B compartments, which are traditionally treated as two weakly related 3D structures and detected with different algorithms. Rather than seeing them as different, our work put them on an equal footing, and we argue that they represent two ends of a continuous spectrum of 3D communities of different sizes.