We propose a new computational method for exploring chromatin structural organization based on Markov State Modelling of Hi-C data represented as an interaction network between genomic loci. A Markov process describes the random walk of a traveling probe in the corresponding energy landscape, mimicking the motion of a biomolecule involved in chromatin function. By studying the metastability of the associated Markov State Model upon annealing, the hierarchical structure of individual chromosomes is observed, and corresponding set of structural partitions is identified at each level of hierarchy. Then, the notion of effective interaction between partitions is derived, delineating the overall topology and architecture of chromosomes. Mapping epigenetic data on the graphs of intra-chromosomal effective interactions helps in understanding how chromosome organization facilitates its function. A sketch of whole-genome interactions obtained from the analysis of 539 partitions from all 23 chromosomes, complemented by distributions of gene expression regulators and epigenetic factors, sheds light on the structure-function relationships in chromatin, delineating chromosomal territories, as well as structural partitions analogous to topologically associating domains and active / passive epigenomic compartments. In addition to the overall genome architecture shown by effective interactions, the affinity between partitions of different chromosomes was analyzed as an indicator of the degree of association between partitions in functionally relevant genomic interactions. The overall static picture of whole-genome interactions obtained with the method presented in this work provides a foundation for chromatin structural reconstruction, for the modelling of chromatin dynamics, and for exploring the regulation of genome function. The algorithms used in this study are implemented in a freely available Python package ChromaWalker (https: //bitbucket. org/ZhenWahTan/chromawalker).
The packing of two meters of DNA in the few-micrometer nucleus results in a structure that performs multiple roles, from forming a structural scaffold of chromatin to facilitating active expression and silencing of genetic material [1,2]. The beginning of interest in the biophysical characterization of chromatin dates to about 50 years ago, spanning from experimental measurements of DNA persistence length [3–5] and thermal stability [4] to pulling individual DNA-protein (DNP) fibrils by convection flows in solution [6], exploring fibril morphology and stability under different media conditions [7], and exposure to ionizing radiation [5]. Before the chromosome conformation capture (3C) [8] era, the classical view of chromatin organization included several successive levels of packing with archetypal structural patterns, ranging from the compaction of nucleosome-bound 10nm fibers [9] with a roughly 200 base-pair periodicity, to the transient 30nm solenoid (hard to detect in vivo) presumably working in the regulation of gene expression [10,11], then to the 30-100kbp loops/domains that are apparently instrumental in shaping large-scale chromatin organization and gene expression [1,12–19]. With the development of the chromosome conformation capture (3C) protocol [8], it has become possible to study chromatin interactions between distant genomic loci. In less than a decade, the original 3C protocol evolved from the analysis of selected pairs of genomic loci to the detection of pairwise interactions between loci and the rest of the genome using chromosome conformation capture on-chip (one-to-all, 4C, [20]), carbon copy (many-to-many, 5C, [21]), and high-throughput 3C (all-to-all, Hi-C, [22]). Finally, improvement of the signal-noise ratio was achieved by performing DNA proximity ligation before nuclear lysis, implemented in in-situ Hi-C [23]. Computational approaches for the analysis of chromatin interaction data developed in recent years can be classified as model-driven or data-driven [24]. Generally, the goal of model-driven studies is to validate physical polymer simulations using Hi-C data. Among them are models of chromatin as a crumpled (fractal) globule [22,25–27], scenarios of loop formation [28,29], analyses of the role of epigenetic factors in driving the chromatin organization [30–34], to name a few. In data-driven studies, on the other hand, experimental Hi-C interaction maps are used for extracting information on statistically significant chromatin interactions, for defining topologically associating domains (TADs) and A/B compartments [35,36], and for the 3D reconstruction of chromatin. Several algorithms have been introduced to study the hierarchical organization of chromatin and its correlation with the distribution of various epigenetic features [37–40], including graph-based approaches for exploring sparse networks of Hi-C interaction peaks, as well as ChIA-PET and HiChIP interaction pairs [41–43]. A recent work by Pancaldi et al. defined chromatin assortativity as a metric for the analysis of correlation between distributions of epigenetic marks and chromatin structure [44]. To date, many methods developed for domain detection [23,45–47] essentially adopt an image segmentation approach aimed at identifying domain regions as a function of short-range interactions along the chromosome, and domain boundary positions are often highly sensitive to the choice of heuristic tuning parameters [48]. Recent network-based methods incorporate effects of long-range interactions in characterizing structural organization [37,39,49] and observe spatial couplings at multiple scales associated with the regulation of gene expression [50]. Spatial reconstructions of chromatin using Hi-C interaction data yield consensus structures [51,52] or ensembles of possible chromosomal conformations [53,54], providing an overall picture of chromatin organization [55]. In this work, we propose a new approach for extracting robust genomic partitions from Hi-C data, seeking to capture the footprints of chromatin structure and organization by considering the entire interaction landscape of this complex system. Specifically, our objectives here are to identify and study structural features of chromatin from Hi-C interaction data and to find a connection between these features and data on epigenetic regulation. Introducing a Markov State Model (MSM) with minimal assumptions and parameters on the chromatin interaction network, we aim to identify structural partitions and interactions between them. By analogy with a biomolecule moving and interacting in condensed chromatin, the MSM allows one to explore chromatin structure using a “probe” randomly walking in the contact energy landscape derived from Hi-C data. Given the multiscale nature of the data-derived contact energy landscape and the metastability of the corresponding MSM, we can identify regions of dense intra- and inter-chromosomal interactions, linkers between these regions, as well as the overall topology of individual chromosomes and the complex structures that chromosomes form by interacting with each other. We found that multiple levels of hierarchy exist in the structure of each chromosome with a layer-by-layer splitting of partitions into subunits with distinct structural and epigenetic features, and presumably, distinct roles. The notion of effective interaction between partitions is introduced and shown to be instrumental in uncovering the hierarchical organization, as well as functional dynamics and epigenetic modulation, of individual chromosomes. Looking at the whole-genome picture, the matrix of effective interactions delineates how chromosomal partitions form a major cluster—with several chromosomes linked by significant inter-chromosomal interactions—as a structural scaffold for genome architecture. The notion of affinity between partitions complements the picture of effective interactions by evaluating the degree of association between partitions, which may contribute to the formation of topologically associated domains, transcription factories and other functional elements, thereby organizing the regulation of genome expression.
A common strategy to study complex network data structures is to combine them with a discrete state Markov process, commonly called a Markov State Model (MSM), with the goal of characterizing hidden network properties [56–58]. MSMs enable one to systematically explore network structure via random walks, where traveling probes form virtual trajectories through the whole network by connecting pairs of nodes. It has been shown that studying the spectral and metastability properties of the network-associated MSM allows one to obtain a reduced description of the underlying complex data. In order to illustrate how MSMs can be used for studying Hi-C data, we introduce here a toy model of a chromosome. Let us consider a linear system characterized by a discrete set of loci S = {1, …, N}, with N = 500. We assume that the number of loci N determines the maximal resolution of this relatively large system. Each locus of the system is associated with an energy Ei, which is linked to the intrinsic stability of the locus i at the given resolution. For the sake of argument, we assume the intrinsic stabilities Ei to follow a hierarchically shaped energy profile (Fig 1A). The energy profile considered here contains 18 wells separated by barriers ranging from 0. 5 to 2 energy units. On the first level of hierarchy, there are two basins separated by a barrier of 2 energy units (black diamonds in Fig 1A), each divided into three sub-basins (indicated as red circles in Fig 1A), which in turn are split into three basins on the third level of hierarchy (black circles in Fig 1A). A traveling probe in such an energy profile is assumed to make two types of moves: sliding between adjacent loci and hopping between non-adjacent ones. We do not make any assumption about the three-dimensional structure of the system and assume a power law contact probability between non-adjacent loci, namely (d0/dij) α, where d0 and dij are the distance between adjacent and any non-adjacent loci respectively, which is equivalent to the genomic distance between loci i and j and such that dij = d0 for adjacent loci. Thus, for each pair of loci i and j, we define the contact energy landscape Eij = (Ei + Ej) /2 − αln d0/dij with α = 1. 5. Assumptions on the power law dependence and the value of exponent α are made on the basis of empirical observations on Hi-C data and polymer models of chromosomes [22]. The contact energy landscape is represented in Fig 1B. To construct the MSM describing the motion of a probe, we define the corresponding Markov generator L for transitions between loci i → j by the Laplacian Lij=e−β (Ej−Ei) /2eβln (d0/dij) α and Lii = −∑j≠iLij, with transition matrix pij = Lij/∑j≠iLij, pii = 0 (∑jpij = 1), flux πij = Lijμj, where steady state probabilities are given by μi=e−βEi/∑je−βEj, and β is an inverse temperature parameter. A network of nodes (loci) and edges (contacts) is obtained from the matrix of fluxes πij, which represents the symmetric probability of contact between a pair of loci. With the set of rules given by the above Markov generator, a probe will tend to explore regions of the network in the neighborhood of the loci that are more stable, i. e. , within an energy well, and will rarely connect loci in different energy wells. This property relates to the “metastability” of the corresponding Markov process. Specifically, in a metastable MSM only a few nodes function as “hubs” of the network, which means that the probe tends to spend most of the time in the neighborhood of these hub nodes, instead of anywhere else. In other words, a probe departing from a generic node in the network is likely to hit the closest hub node in the set hub-nodes M. Additionally, the probability for a probe departing from a hub-node in M to return to itself is larger than that for the probe to reach another hub-node in M. As a result, nodes in the neighborhood of hubs tend to cluster together in a modular manner. This is a condition that allows one to find a reduced size MSM that approximates the original Markov process associated with the initial network. One can quantify how well the probe motion satisfies this condition by defining a metastability index ρM. A metastability index is the ratio of two probabilities (see Eq 6 in Methods for a precise definition): the probability Pout for a probe to connect two different hubs in the set M (as small as possible) over the probability Pin for a walker to hit any hub in the set M irrespective of the starting point (as large as possible) [59]. In a metastable MSM the metastability index is expected to be a small number (ρM=Pout/Pin<1) characteristic of the hub set M. To understand how metastability works, it is instructive to consider the large changes in kinetic properties of the MSM upon increase of the inverse temperature parameter β (annealing condition). These changes are clearly illustrated by the mean first passage time MFPT τij, which is the average time (number of steps) a probe takes to connect the pair of states i and j (where states represent loci of the toy chromosome, see Eq 4 in Methods). Fig 1C shows the MFPT matrices in the case of low β = 1 and high β = 10, respectively. A clear separation of time scales emerges upon increasing 𝛽, as reflected in the partitioning of the MFPT matrices. The nested squares emerging in the MFPT matrix (Fig 1C) at high β identify pairs of states/loci (i, j) with comparable values of the MFPT τij, which is a result of the hierarchically shaped energy profile. As β is increased, the emerging separation of time scales in the MFPT matrix is the result of the dominant barrier that separates a given pair of loci in the energy landscape (see Fig 1A): each of the separated regions contains one or more hubs that cause probes to stay within its vicinity. As a result, the effect of high β on the MFPTs of the MSM elucidates how dominant interactions in the system can be captured using just a subset of loci, the hub set M. To quantitatively identify the hub set, an optimization procedure is performed in order to find the sets M that minimize the metastability index ρM (see details in Materials and Methods) as a function of increasing β. Fig 1D shows the optimized profile of the index ρM as a function of the hub set sizes, and at different values of β. All the profiles of ρM clearly show three minima corresponding to the hub sets M (2), M (6), and M (18) (of sizes 2,6 and 18, respectively), which correctly identify locations of the energy wells in the hierarchically shaped energy landscape in Fig 1A. The hub sets obtained by optimizing the index ρM are suitable as cores of partitions, which characterize the coarse-grained state space of an approximated MSM. Fig 1E (top) depicts the toy network associated with the contact energy landscape shown in Fig 1B. Nodes are colored according to the partitions constructed around nodes in the hub set M (6). A reduced network corresponding to the hub set M (6) is also shown in Fig 1E (bottom). The nodes in this network are defined as soft partitions of the initial set of loci S, whereas the links characterize the “effective interactions” between nodes with values Fab = ∑i∈Sqa (i) πib (see Eq 11 in Materials and Methods). The quantity qa (i) is a committor probability, which is the probability for a probe departing from a locus i to hit the locus a∈M before any other locus in the hub set M (see Eq 8 in Materials and Methods). Using the intuition acquired with the help of this toy model, we describe in the following section how a MSM can be constructed from the Hi-C dataset of a single chromosome and how metastability analysis can be performed in order to infer chromosomal architecture and effective interactions between partitions. We now consider the random walk through the interaction network of a single chromosome, using the example of Hi-C data on human chromosome 17 in the human B lymphoblastoid cell line GM12878 at 50kbp resolution [23] and describing it via a Markov process. To do that, we start from the number of times fij a pair of genomic loci i and j is found in a contact. After applying a Gaussian smoothing filter on the raw data (see Hi-C data preprocessing in Materials and Methods), a pairwise contact energy Eij = −lnfij is defined for each pair of genomic loci. With this interpretation, the larger the contact frequency the more stable (lower contact energy) pair of genomic loci is involved. The representation of this two-dimensional contact energy landscape is shown in Fig 2A. A probe moving in such a landscape is expected to spend most of the time in pairs characterized with low contact energy and rarely connecting across high contact energy pairs. In the toy model presented in the previous section, a pairwise contact energy landscape (Fig 1B) was constructed from the one-dimensional energy landscape (Fig 1A). Here, we use a reverse logic and consider the one-dimensional projection (Fig 2B) of the two-dimensional contact energy landscape (Fig 2A). To do that we define the contact energy of a genomic locus i as Ei = −lnfi, where fi = ∑jfij is the total number of times a genomic locus i is found in any contact, hence loci involved in more contacts are more stable as they exhibit lower contact energy. Fig 2B shows the 1D projection of the pairwise contact energy landscape (for both raw and Gaussian-smoothed data), which presents multiple features—minima, maxima, and barriers—characterizing the architecture of the chromosome. Here, we briefly describe the metastability analysis applied to chromosome 17 (steps 1–4) and consider whole-genome interactions (step 5) using a coarse-grained approximation. Step 1. In order to explore chromosomal architecture, the MSM describing the motion of a probe in the contact energy landscape Eij is implemented by introducing the Maxwell-Boltzmann probability πij (β) =e−βEij/Z (β), where Z (β) =∑ (i, j) e−βEij is the partition function, β is the inverse temperature parameter, and πij is the symmetric flux of probes between pairs of loci. Using a 50kbp resolution for the Hi-C dataset, a total of N = 1625 genomic loci comprise the state space S of the MSM for chromosome 17. The transition matrix associated with the MSM is defined as pij = πij/μi, where μi (β) = ∑jπij (β) is the Boltzmann weighted probability (steady state probability distribution) of observing locus i involved in any contact. The effect of annealing (increasing the inverse temperature β) on the kinetics of a random walker is clearly reflected in the MFPT matrices (Fig 2C), obtained at low and high β, respectively. While at low β (β = 1) MFPTs show no partitioning, a separation of time scales becomes evident at high β (β = 9). Indeed, the 1D projection of the pairwise pseudo-energy landscape (Fig 2B) shows that, apart from the centromere that naturally separates the two chromosome arms, the highest barrier in the 1D projection is about 1. 5 in β−1 units (see Fig 2B). Therefore, partitioning of MFPTs scales is observed only for significantly higher values of β. Step 2. Optimization of the metastability index ρM (see details in Materials and Methods) over the hub set M of different sizes was performed as a function of the inverse temperature parameter β, revealing the levels of structural hierarchy of chromosome 17. The ρM profiles upon increasing β (Fig 2D) converge towards five minima, which correspond to the hub sets M (2), M (5), M (8), M (12), and M (27), of sizes 2,5, 8,12, and 27, respectively. The M (2) hub set is not considered as it trivially identifies the chromosome arms separated by the centromere. It should be noted that the locations of the obtained hub sets correspond to the locations of the multiple wells present in the projected contact energy landscape, as shown in Fig 2B. Step 3. Given the hub sets obtained at different levels of structural hierarchy, one can identify chromosomal partitions, namely regions of the chromosome compacted around corresponding hubs and, at the same time, separated from one another. Soft partitions are defined around corresponding hubs using the committor probability qa (i) [59], which in this case is interpreted as the probability for a locus i to belong to the partition defined by the hub a∈M. To identify physical partitions of the chromosome in relation to other chromosomes, a coarse-grained description is adopted here by considering hard partitions. In this case, a step function θA (i) characterizes whether a locus i belongs to a partition A, specifically θA (i) = 1 if i ∈ A, θA (i) = 0 otherwise, and ∑AθA (i) = 1 for any locus i (see Eq 9 in Materials and Methods). Fig 2E illustrates the partitioning of the network for human chromosome 17 that is obtained from the hub set M (12). Step 4. To complete the description of chromosome structure, one needs also to characterize the strength of interactions between the partitions obtained at different levels of hierarchy. As in the example illustrated in the toy model, we consider the effective interaction between two soft partitions located around the hub loci a and b of a chromosome as the mean contact energy acting between them, which corresponds to the weighted flux connecting loci a and b via the committor probability qa (i), namely Fab = ∑i∈cqa (i) πib (see Eq 11 in Methods). Step 5. In the context of whole-genome interactions, a coarse-grained description is adopted (see Step 3) for estimating the mean contact energy between pairs of partitions in the 23 chromosomes: FAB = ∑i∈gθA (i) ∑j∈gπijθB (j), where θA (i) and θB (i) are step functions and πij is the flux of probes between corresponding loci (see Materials and Methods for details). Fig 3A–3C show the partitioning of chromosome 17 at three levels of hierarchy with corresponding effective interaction matrices (Fig 3D–3F), and the band representation of partitions at all three levels (M (5), M (12), M (27); Fig 3G). The major partition boundaries that emerge at the first level of hierarchy persist through the other levels (Fig 3G; similar for all chromosomes, see S1 Fig) and show a qualitative agreement with the borders of euchromatic and heterochromatic bands obtained from Giemsa staining (Fig 3G and S1 Fig). Unfortunately, as Giemsa staining is a very basic and crude cytological method for identifying densely-packed (heterochromatic, dark stain) and low-density (euchromatic, light stain) genomic regions, it is not possible to perform an accurate quantitative analysis on staining bands [61,62]. At the lowest level of hierarchy (Fig 3A), we observed the bulk topology of the chromosome where two chromosomal arms are brought together via strong interactions between partitions 1,4, and 5. Most partitions at the lowest levels of hierarchy are found to contain both euchromatic and heterochromatic bands, and they have highly distinct structural and/or functional characteristics. At the second level of hierarchy (Fig 3B), partitions 1. 2,4. 1,4. 2,5. 1, and 5. 4 form several non-adjacent contacts, working as hubs responsible for most of the network structure. The third level of hierarchy (Fig 3C) yields further details of chromosomal architecture: the p-arm is loosely connected and is weakly centered on 1. 2. 2 and 1. 2. 4, while the q-arm is densely connected by multiple hubs (4. 1. 1,4. 2. 2,4. 3. 1,5. 1. 2, and 5. 4. 1). At this level many partitions are homogeneous, either eu- or heterochromatic, interacting more strongly with partitions with similar packing densities, resembling the phenomenology of the so-called A/B (active/inactive) chromatin compartments [22]. For instance, partition 1. 2 is split into mostly euchromatic (1. 2. 2,1. 2. 3) and heterochromatic (1. 2. 1,1. 2. 4) partitions, while partition 4. 1 is split into predominantly euchromatic (4. 1. 1,4. 1. 2,4. 1. 3) and heterochromatic (4. 1. 4,4. 1. 5) ones. Partition 4. 1. 1 is the largest among these, forming significant interactions with the p-arm through partition 1. 2. 2. Another noticeable interaction between chromosomal arms occurs via the partition 4. 3. 1, which links heterochromatic partitions 5. 1. 1–2 and 1. 2. 4. Interestingly, the mostly euchromatic partition 1. 2. 2 is responsible for many non-adjacent contacts with the q-arm, whereas heterochromatic 1. 2. 4 forms non-adjacent contacts only with 4. 3. 1 and 3. 1. 1. With these observations, one may conclude that heterochromatic partition 1. 2. 4 acts as a structural foundation that link the mostly euchromatic partitions 1. 2. 2,1. 2. 3,2. 1. 1,2. 1. 2, and 3. 1. 1. To investigate how the hierarchical organization of chromosomes facilitates their function, we first analyzed the average density of various epigenetic factors in partitions (Fig 4 and S3 Fig), using chromosome 17 as an illustration for this analysis and operating at the third level of structural hierarchy. Fig 4A, in which node sizes depict partition sizes, shows that heterochromatic partitions 5. 1. 1 and 5. 1. 2 apparently form a structural foundation of chromosome 17 architecture, linking the p- and q-arms through the large mixed partition 4. 3. 1 and the heterochromatic partition 1. 2. 4. Next, we consider two transcription factors commonly associated with chromatin structure studies, namely CTCF (transcriptional repressor, Fig 4B) and RAD21 (cohesin, S3H Fig). The CTCF graph (Fig 4B) shows that the heterochromatic partitions (1. 2. 4,5. 1. 1,5. 1. 3, and 5. 2. 2) and the pericentromeric partition 3. 1. 1 have the lowest CTCF levels, while the highest CTCF levels were found on 4. 2. 2,4. 2. 1, and 5. 4. 3. The euchromatic or mostly euchromatic partitions 1. 1. 1,1. 2. 2,2. 1. 2,2. 2. 1,4. 2. 3,5. 4. 1, and 5. 4. 2 show average levels of CTCF in the overall eight-fold variation in the density of this transcription factor across partitions. Among the hub partitions, namely those that form extensive non-adjacent contacts, only 4. 2. 2 shows high CTCF levels. Unlike CTCF, RAD21 (a component of cohesin) exhibits only a two-fold variation in densities across partitions at this level of hierarchy. The correlation between CTCF and cohesin binding sites has been noted previously [63,64], and indeed the distribution of RAD21 (S3H Fig) chiefly follows the same general trends as that of CTCF. The strongest among the few exceptions are the increased density of RAD21 in 4. 2. 3 and decreased density in 5. 4. 3. Turning to histone modifications, we note that H3K9ac (Fig 4C) and H3K9me3 (S3A Fig) are associated with activation and silencing of transcription in corresponding promoter regions and, therefore, are expected to show opposite density trends. Indeed, densely packed heterochromatic partitions (1. 2. 4,5. 1. 1,5. 1. 3) and pericentromeric 3. 1. 1 show very low levels of the activating H3K9ac histone modification, while the silencing H3K9me3 modification shows increased density in these partitions (highest in the case of 3. 1. 1). At the same time, euchromatic and mostly euchromatic partitions 1. 1. 1,1. 2. 2,2. 1. 2,4. 1. 2,4. 2. 2,5. 4. 2, and 5. 4. 3 show an increased density of both epigenetic factors, with some slight variations. The opposing trends are observed in heterochromatic partitions for the activating H3K27ac (decreased density, S3B Fig) and inhibiting H3K27me3 (increased density, S3C Fig) modifications, with the most pronounced effects being on 1. 2. 1,4. 1. 4,4. 1. 5,5. 1. 1, and 5. 1. 3. Distributions of the H3K4me1 (S3D Fig) and H3K4me3 (S3E Fig) modifications—both activators—show higher densities in most euchromatic partitions, and in few heterochromatic ones—1. 2. 1,4. 2. 1, and 5. 1. 2. Interestingly, the heterochromatic partitions 1. 2. 1,4. 2. 1, and 5. 1. 2 are enriched in all activating histone modifications considered here (H3K4me1, H3K4me3, H3K9ac, H3K27ac), and, at the same time, are depleted in the inhibiting modifications H3K9me3, H3K27me3. These trends suggest that the above partitions may contain facultative heterochromatin that switches between active and repressed states. Overall, the DNA accessibility graph, indicating the DNase-Seq signal (Fig 4D), shows that most of the euchromatic partitions (1. 1. 1,1. 2. 2,2. 1. 2,4. 2. 2,5. 4. 2, and 5. 4. 3) are rather open and accessible for contacts or interactions. Increased accessibility observed for partitions 4. 2. 1 and 1. 2. 1 is consistent with the conclusion that these partitions may contain facultative heterochromatin, which was inferred from the distribution of activating and inhibiting histone modifications. The partition 5. 1. 2, on the contrary, is less accessible, suggesting that it contributes mostly to the structure formation. Finally, the distributions of RNA polymerases II and III (S3F and S3G Fig) complement the picture of the potential functional involvement of different partitions in chromosome 17. RNA polymerase II (POL2), crucial component of mRNA synthesis, is distributed quite evenly in both euchromatic and heterochromatic partitions (except the high level in 2. 2. 1). The synthesis of tRNA, 5S rRNA, and small RNAs through the action of RNA Polymerase III (POL3) is distributed in a more specific way across different partitions. The POL3 signal is high in euchromatic partitions 1. 2. 2,4. 1. 2,4. 2. 3, in mixed 4. 2. 1 and 5. 3. 1, as well as in some heterochromatic ones (4. 1. 4,4. 1. 5,5. 1. 2,5. 1. 3, and 5. 2. 2). Peculiarities in distributions of epigenetic factors, DNA accessibility, and RNA polymerases revealed in the analysis of individual chromosomes should be further considered in the framework of whole-genome organization, exploring the interplay between intra- and inter-chromosomal interactions in the regulation of gene expression. To this end, we moved from single-chromosome analysis to studying the whole-genome effective interaction matrix. Given that chromosomes are spatially segregated into chromosomal territories (CTs), one can approximate the whole-genome organization by merging single-chromosome partitioning schemes at appropriate levels. Using a selected representative level from each chromosome (see Materials and Methods: Chromosome partitioning), we formed a whole-genome description with 539 partitions, with an average partition size of about 5Mbp (S2 Table). The matrix of effective interactions between chromosomal partitions (Fig 5) provides a general view of the overall physical interactions in chromatin. It shows that chromosome 1 and small chromosomes (14–20 and 22) massively interact with others, while chromosomes 4,5, 9,21, and X appear to be relatively isolated from the rest of the genome (Figs 5 and 6). Several partitions form consistently stronger intra- and inter-chromosomal interactions with other partitions. We classified interaction strengths into 5 layers with equally-spaced threshold values: the scaffold layer is the strongest, followed by layers 1,2, etc. (see also Materials and Methods for the definition of the interaction strength at different layers). Fig 6A shows the major cluster in the whole-genome partition set: partitions from different chromosomes form tight sub-clusters highlighted by color and marked by chromosome labels. All displayed partitions are linked by the two strongest layers of interactions (scaffold interactions are represented by black edges, and layer 1 by grey edges). It is easy to see that most of the intra-chromosomal contacts and some inter-chromosomal interactions are established on the scaffold layer, giving rise to a structural foundation for genome-wide architecture (Fig 6). Specifically, chromosomes 1,14,16,17,19,20, and 22 are densely interconnected, while other chromosomes in the major cluster are linked to them via only a few interactions. Partitions 1–2. 1. 2,14–3. 4, and 22–6, for example, act as contact hubs between these massively interacting chromosomes and others. On the other hand, partitions such as 3–2. 2. 2,8–6. 2. 1,10–3. 4. 1, connect less-strongly interacting chromosomes to the strongly interacting ones (see S3 Table for interaction strength layers for these interactions between chromosomes). Notably, chromosomes 1 and 2, the two largest ones (about 250Mbp each), behave differently in the context of the whole-genome interactions. While chromosome 1 serves as a hub in the interactions between the highly- and less-interacting chromosomes, chromosome 2 does not show many interactions with other chromosomes (Fig 6A). Turning to functional regulation, most of the partitions involved in significant inter-chromosomal interactions exhibit higher densities of several epigenetic factors, such as CTCF (Fig 6B), H3K9ac (S5A Fig), and DNase accessibility (S5B Fig). These partitions may participate in the formation of active epigenetic compartments facilitated by the structural role of CTCF [65]. Active processing of genomic information taking place in these structures is regulated by the activating histone modifications (H3K9ac) and transcriptional repressors (CTCF). The opposite trend is observed for partition 14–3. 4, which is coupled with a higher density of the silencing H3K9me3 histone modification. Therefore, partition 14–3. 4 and its interactions with partitions in other chromosomes, for instance 10–3. 4. 1, with low activating factor densities may indicate the formation of dense structural heterochromatin and/or silencing facilitated by Polycomb bodies [1,66]. To evaluate how the distribution of epigenetic signals may be associated with interaction between partitions, we calculated correlations between effective interaction strengths and the expected enrichment of factor densities across partition pairs that are mostly euchromatic (EC) or heterochromatic (HC). The enrichment of factor densities is estimated here as the product of factor densities per partition (S12 Fig). To obtain the strongest signals, we limited our consideration to interactions between partitions that are dominated by either eu- or heterochromatin (see legend for S12 Fig for the definition of EC and HC partitions): EC-EC pairs (S12A Fig), HC-HC pairs (S12B Fig), and EC-HC pairs (S12C Fig). Despite the relatively weak correlations, general trends appear to be quite clear, with the strongest ones seen between euchromatic (EC-EC) partitions (S12A Fig). Transcription factors CTCF and RAD21 are always positively correlated, as well as POL2 in EC-EC (S12A Fig) and EC-HC (S12C Fig) pairs, whereas POL3 shows no correlation. The positive correlation for CTCF and RAD21 with effective interaction strength agrees with current literature on the role of CTCF and cohesin in mediating chromatin structure through looping interactions [65,67,68]. Stronger interactions between EC partitions appear to be linked to higher transcriptional activity, as suggested by the positive correlation with active histone modifications and POL2. Absence of correlation for HC-HC pairs in the case of POL2 can be related to the fact that transcriptional activity is suppressed in heterochromatin. Potential active involvement of interacting euchromatic partitions in the formation of transcription factories is corroborated by the most pronounced correlation observed for DNA accessibility in pairs of euchromatic partitions (S12A Fig). Activating histone modifications, except for H3K4me3, show positive correlations in all types of interacting partition pairs. Interestingly, silencing histone modifications appear also to be weakly correlated with effective interactions between partitions. The original partitioning analysis was performed on the GM12878_primary (B lymphoblastoid) Hi-C dataset by Rao et al. [23] (GEO accession GSE63525). We also applied our model to four other datasets: GM12878_replicate (a biological replicate of GM12878_primary dataset), IMR90 (lung fibroblast), HUVEC (umbilical vein endothelial cells), and HMEC (mammary epithelial cells). Our goal in this analysis was two-fold: (i) to benchmark robustness and reproducibility of the method using the replicate dataset; (ii) to examine the sensitivity of the method in detecting alterations in chromatin organization in different cell lines, associated with corresponding genome functional states and gene expression levels. S13 Fig shows side-by-side comparisons of the partitioning schemes for GM12878_primary and the other datasets, and S6 Table shows some indicative statistics comparing the results from each case. First, we observed a high consistency between the biological replicates of GM12878: S13A Fig shows that the partitioning was highly consistent between the two sets of Hi-C data, with partition boundaries being identical in most cases, resulting in the high Rescaled Mutual Information (RMI) of 0. 70 (see Chromosome partitioning in Materials and Methods for definition of RMI). The composition of the major cluster was also largely identical. Comparing the results from other cell lines, we observed significant differences: IMR90, HUVEC and HMEC cells each had significantly shifted partition boundaries compared to GM12878_primary, leading to lower RMI values of 0. 39 to 0. 48. The major-cluster structures in these cell lines are also significantly different (see S16 Fig for IMR90 and HUVEC), especially that of HMEC, where no strong inter-chromosomal interactions were observed between partitions, and the chromosomes remained isolated in the whole-genome network. Notably, in both IMR90 and HUVEC, a large partition on chromosome 9 forms extensive inter-chromosomal interactions: the overlapping region (chr9: 1268000000–1412500000) contains two genes (OLFM1 and MVB12B) with the RNA-expression profiles different from that of GM12878. The MVB12B (a component of endocytic protein system [69]) gene is activated in both IMR90 (lung fibroblast) and HUVEC (umbilical vein endothelial cells) cell lines, and OLFM1 (lung cancer marker [70]) in IMR90, while both genes are silenced in GM12878. These preliminary observations call for future in-depth investigation of the structural basis, functional mechanisms, and specifics of epigenetic regulation behind the observed differences between cell types. While effective interactions between partitions characterize the overall architecture of genome organization, it may not fully discriminate functionally relevant interactions between chromosomes and their parts. Indeed, most partitions are presumably in constant motion within the nucleus, and as Hi-C experiments are typically conducted on unsynchronized cell populations, effective interactions capture the average contact probability arising from both random diffusion and specific transient interactions. Therefore, in addition to effective interactions, the affinity between partitions was also calculated, which reflects how the observed interaction frequency differs from the expected frequency (from random diffusion), because of possible associations between partitions. Defined as the ratio between observed and expected contact probabilities between pairs of partitions (see Eq 15 in Methods), the affinity is indicative of the degree of association between partitions, and high affinity values may serve as a manifestation of biologically-relevant contacts. Fig 7 contains the whole-genome matrix of pairwise affinities (blue: high affinity, white: low affinity) between corresponding partitions. Like the observations in the whole-genome effective interaction matrix (Fig 5), the largest chromosomes 1 and 2 exhibit different behavior, with chromosome 1 containing partitions with high affinity to those in several other chromosomes (especially with chromosomes 14–22) and chromosome 2 generally showing low affinity to partitions in other chromosomes. Smaller chromosomes 14–22 form more, presumably functional, contacts with each other, compared to other chromosomes. At the same time, the number of partition pairs with high affinity is much lower than number of pairs with significant effective interactions (compare Fig 5 and Fig 7). In total, we observed 687 high-affinity pairs (S4 Table), which are seemingly crucial for whole-genome structural organization and function. Interestingly, several large partition pairs (>2Mb) with high effective interactions and affinity are located in the telomeric regions of corresponding chromosomes (yellow cells in S5A Table), having moderately high densities of epigenetic/transcription factors and increased DNA accessibility (S5A Table). Two other groups of partitions with high affinities are characterized by smaller partition sizes and highly elevated concentrations of various transcription factors and epigenetic modifications (S5B Table): (i) pericentromeric partitions (red cells in Table) show high concentrations of activating (H3K4me3) and silencing (H3K27me3 and H3K9me3) histone modifications and high levels of POL3 and RAD21; (ii) telomeric partitions (yellow cells in Table) show strongly increased concentrations of all activating histone modifications, POL2, and CTCF, as well as high DNA accessibility. This separation between types of activating histone modifications, transcription factors, and DNA accessibility in centromeric and telomeric regions signals a specificity of functional interactions between partitions with high affinities to each other. Examples of partitions involved in high-affinity interactions and characterized by the over-representation of different epigenetic factors and modifications are collected in Fig 8 and S6 Fig, where high-affinity clusters of partitions enriched in these epigenetic marks are plotted. A comparison of the inter-chromosomal interactions in the major cluster of effective interactions (Fig 6) with interactions in affinity clusters (Fig 8 and S6 Fig) highlights several relatively-small partitions, e. g. 9–5. 6. 3,9–5. 6. 5, and 11–5. 1. 3, that act as junctures between different chromosomes. These partitions yield increased density of CTCF along with other juncture-partitions (3–2. 2. 2,6–1. 4. 3,8–6. 2. 1 to name a few), pointing to the potential importance of these partitions in whole-genome structural organization. This inference is further supported by multiple interactions detected for partitions 8–6. 2. 1,9–5. 6. 5, and 11–5. 1. 3 in the CTCF affinity graph (S6D Fig). Additionally, H3K9ac (Fig 8A) and H3K27ac (Fig 8C) affinity graphs hint at the functional importance of some of these partitions: the central part of the H3K9ac graph is formed by partitions 9–5. 6. 3,9–5. 6. 5, and 11–5. 1. 3, while 9–5. 6. 3 and 11–5. 1. 3 are also present in the H3K27ac graph. Focusing on individual epigenetic factors, the activating H3K4me3 mark links more partitions than the activating H3K4me1 histone modification. The silencing H3K9me3 histone modification functionally links many centromeric partitions, whereas the activating H3K9ac modification works in both mostly euchromatic and mixed euchromatic/weakly-heterochromatic non-centromeric regions. Similarly, the activating H3K27ac modification affects mostly non-centromeric partitions, unlike the very active silencing H3K27me3, for example, in partitions 1–4. 11. 1,2–3. 1. 1, and 10–2. 1. 1 (Fig 8). These partitions are also characterized by the high levels of POL3 (S6B Fig) and RAD21 (cohesin, S6C Fig), whereas the insulator CTCF links several euchromatic partitions across different chromosomes (S6D Fig). It is evident that centromeric partitions 1–4. 11. 1,2–3. 1. 1, and 10–2. 1. 1 are enriched with almost all regulatory factors (see Fig 8 and S6 Fig), yielding high affinities to other partitions and pointing to important functional interactions and intense regulation taking place in these partitions. Interestingly, while activating histone marks (Fig 8A, 8C, 8E and 8F) are dominant in several euchromatic partitions, these marks are also present in partitions containing large sections of heterochromatin and centromeres, which are commonly associated with dense packing and transcriptional repression. Similarly, silencing histone marks (Fig 8B and 8D) are dominant not only in heterochromatic and centromeric partitions, but also in some partitions that are mostly euchromatic. Furthermore, dominating regions for the transcription factors CTCF and cohesin (S6C and S6D Fig) appear to have significant overlap with activating and silencing histone marks, respectively. These overlaps show the complexity of functional interactions in chromatin, even at the coarse-grained level of partitions: opposing factors are found acting in the same regions, allowing for switching between transcriptional states in response to other biochemical cues.
We proposed here a computational framework for exploring chromatin organization based on Markov State Modelling of chromatin interactions. Given the multilevel hierarchical packing of chromatin, we introduced a reduced description of the complex network of chromatin interactions and its organization via interactions between structural units at different levels of hierarchy. By interpreting Hi-C data as a pairwise contact energy landscape, a Markov State Model approach was used to explore the chromatin interaction network through the random walk of a probe. While steady-state distributions obtained from the Markov process of randomly-moving molecules can serve as a measure of the chromatin accessibility for epigenetic factors [71], taken alone they describe neither the genome architecture, nor structural and functional interactions between genome partitions and regulatory factors. In this work, analysis of the Markov State Model under thermal annealing shows the key role played by the ruggedness of the contact energy landscape in shaping chromosome structural organization. Specifically, metastability analysis of the Markov State Model associated with the chromatin interaction network allowed us to identify levels of structural hierarchy and to observe structural units—partitions of different scales. These structural partitions serve as a coarse-grained description of chromosomes, which form the basis for introducing the notion of intra- and inter-chromosomal network of effective interactions. The analysis of effective interaction networks across levels of hierarchy in individual chromosomes shows that chromosomes adopt highly varied topologies. While the lower levels reveal an overall architecture of the folded chromosome, the higher levels can provide structural details in relation to functional organization and regulation of gene expression. Biological insight on the structural organization of chromosomes can be obtained with our method by considering peculiarities in the distributions of transcription and epigenetic factors in eu- and heterochromatic partitions in relation to interactions between them. Partitions at the highest levels of hierarchy may be seen as analogous to TADs, or to the so-called A/B (active/inactive) epigenomic compartments [72]. In the future, with the development of common standards and benchmarks by the community, it would be important to compare results and insight obtained from various genomic segmentation approaches. In this work, however, we based our analysis on the obtained sets of partitions, showing how studying distributions of activating and silencing histone modifications in these partitions can help to understand the role of structural organization of chromosomes in the regulation of gene expression. Shifting our focus from structural analysis of individual chromosomes and the functional involvement of partitions to whole-genome architecture, we considered the set of 539 partitions obtained at high levels of hierarchy in corresponding chromosomes with effective interactions between them, which were obtained by adopting a fast “mean field” approximation. In the context of the “partition space”, an analysis of genome-wide effective interactions provides a blueprint of inter-chromosomal contacts, showing that despite the strong crowding of partitions in chromosome territories, most chromosomes are significantly connected with each other, giving rise to a bulky cluster in the core of the effective interaction network. The strongest interactions were observed between chromosomes 14–22, which are characterized by small chromosome sizes. Heterochromatic partitions are apparently mostly involved in the formation of chromosomal territories, interacting within the corresponding chromosomes and providing structural integrity, and showing only low levels of activating factors. On the other hand, most of the inter-chromosomal juncture partitions, while relatively small in size, are enriched with CTCF, H3K9ac, and DNase-Seq, which may lead one to conclude that these partitions are involved in the formation of inter-chromosomal active epigenomic compartments [72]. Correlations of effective interactions between partitions with distributions of epigenetic factors in these partitions show that: (i) most active regulation apparently takes place in pairs of interacting euchromatic partitions; (ii) DNA accessibility, CTCF and activating histone modifications H3K4me1, H3K9ac, and H3K27ac are major potential contributors in the regulation of genome function. Additional biological insight was obtained by determining partitions that may form transient function-related contacts, thereby triggering alternate chromatin states. To this end, the affinity measure was introduced here to evaluate the level of association between partitions, so as to identify partitions with functionally-related interactions. Irrespective of the effective interaction strength, high affinity between partitions point to the mutual functional involvement of corresponding partitions. Since different factors are likely to play dominant roles in different genomic regions, our affinity analysis is complementary to the concept of chromatin assortativity introduced by Pancaldi et al. [44], which may identify epigenetic factors associated with multiple high-affinity communities across the whole genome. There are different challenges in extending the original analyses of the Hi-C data to exploring the structure-function relationships in the genome. Several previous studies using a hierarchical clustering approach for the analysis of Hi-C data are based on the a priori assumption of the existence of structural hierarchy in chromatin [38–40]. While the work by Boulos et al. is free from such an assumption [37], it employs a tunable scale parameter in establishing the hierarchy. Our approach is based on an energy landscape representation of the chromatin interaction network, which is formalized and explored via a Markov State Model. Metastability analysis of the Markov State Model allows one to detect natural levels of hierarchy in chromatin structure. The novelty of our approach does not free it, however, from certain limitations. For example, the “mean-field” approximation for coarse-graining of whole-genome interactions, rooted in the computational challenge in calculating committor probabilities on very large networks, is currently a necessary step for processing massive whole-genome datasets. On the other hand, because of its rapid computational time, this approximation allows us to explore higher resolution Hi-C datasets and to obtain partitions of smaller sizes and, ultimately, uncovering higher levels of hierarchy. Further, as the rugged morphology of the energy landscape is key to determining chromosomal partitions, the effect of noise associated with the Hi-C data is of critical importance. For example, all contact frequencies and derived quantities are affected by the characteristic noise, since most Hi-C experiments are conducted on unsynchronized cell populations. While a possible solution would be to study single-cell Hi-C [73] data, which shows great promise in capturing differences between transient states in chromatin organization, the current protocol yields too few interaction pairs for a meaningful analysis of the interaction network. Although specialized variants of the Hi-C protocol, such as capture Hi-C (cHi-C) [74], do not provide a full view of the physical organization of chromatin, they can nonetheless be useful in targeting specific subsets of genomic loci, such as promoters. The biological implications of our analysis, particularly on the relationship between factor enrichment and effective interaction strength and affinity, may also be strengthened by incorporating additional experimental approaches, such as ChIA-PET [75] and HiChIP [76], which identify interacting genomic elements that are concurrently associated with specific binding proteins. Genome architecture mapping (GAM), a newly devised experimental protocol that determines the frequency at which genomic loci lie on the same spatial plane by sequencing fragments isolated in cryosections of the nucleus [77], can be a great source of constraints for future 3D whole-genome reconstruction. To conclude, there is no doubt that scientific interest in chromatin structure will continue to drive the development of a variety of specialized experiments and computational approaches in the field of 3D genomics. The method presented here, aimed at detecting and characterizing the hierarchical organization of chromatin, is a step towards unravelling causal relationships in chromatin structure and dynamics of function-related transient molecular phenotypes. The great potential of new experimental data combined with constant methodological improvement are critical in the quest for a more detailed understanding of chromatin architecture, 3D reconstruction, dynamics, and epigenetic regulation.
No human or animal subjects and/or tissue were used in the work. Ethics rules of the Bioinformatics Institute, A*STAR were followed during the work on the project and preparation of the paper. In the following, a Markov jump process is introduced to describe a random walk in the chromatin interaction network, where a probe connects pairs of interacting genomic loci, which represent the states of the Markov State Model (MSM). We first focus on a single chromosome c and denote the corresponding matrix element of Hi-C counts fij for a pair of loci (i, j). We define a pairwise interaction pseudo-energy Eij = −lnfij, which characterizes a strength of interaction between a pair of loci (i, j): the higher the observed counts the more stable the corresponding interactions (lower pseudo-energy) are. Therefore, a Maxwell-Boltzmann probability distribution of counts is defined as πij (β) =1Z (β) exp (−βEij) (1) where Z (β) = ∑ (i, j) ∈c exp (−BEij) is the partition function and β is the thermal parameter (inverse of a temperature). Eq 1 evaluates the joint interaction probability of a pair of loci (i, j), which can be modulated by the thermal parameter β. For low values of β (high temperature), the interaction energies tend to contribute equally in the exponential, whereas for high values β (annealing, low temperature) only the highly interacting loci contribute significantly in the exponential. The transition probability associated with the Markov jump process is defined by the following conditional probability pij=πijμi (2) where μi=∑j∈cπij (3) is the probability for the genomic locus i to form any interaction in chromosome c. Additionally, μi is the steady state distribution of the transition matrix pij. The transition probability matrix in Eq 2 uniquely identifies a discrete time Markov jump process that governs the trajectories of a random walker across the state space. A random walker is interpreted as a probe particle traveling between genomic loci, for instance a protein such as a transcription factor. Accordingly, the steady state distribution μi can be interpreted as a distribution of probes in a locus i, whereas the distribution πijc in Eq 1 describes an undirected flux of probes connecting the loci i and j. The kinetic distance between pairs of loci is the mean first passage time (MFPT), the mean number of discrete steps τij between two different genomic loci (i, j) in chromosome c, which is obtained by solving the system of equations [78] τij=pij+∑k≠i, jpik (1+τkj) (4) If the departure and arrival states coincide, the MFPT is called mean recurrence time MRT τi, which gives the mean time for a walker to return to its initial state i. The MRT is obtained from the MFPTs in Eq 4 via the formula τi=pii+∑k≠ipik (1+τki) (5) S7 Fig shows the MFPT matrices (with the MRTs in diagonal) for chromosomes 1,17, and 20 at β = 1 (left) compared to annealing condition at high β (right). A separation of time scales emerges upon increasing the β parameter, which is reflected in the partitioning of the MFPT matrices. The squares depicted in the annealed MFPT matrices identify sets of pairs (i, j) with a similar τij, which are the results of the partitioning in the network of interactions. S8 Fig shows the MFPT matrices for all chromosomes at corresponding high β. For each chromosome, the value of β was chosen as large as possible to observe the fine partitioning structure in the chromosome, while at the same time avoiding singular values in the calculations of the MFPTs and MRTs. The values of β used are listed in S1 Table. The optimization of the metastability index upon annealing conditions allows one to obtain an optimal hub set M. It has been shown that in a MSM with a large state space, an optimal hub set and its corresponding partitions can be used for reducing the state space and obtaining a smaller MSM [59]. In the context of our model of chromosome interactions, we introduce a scheme for coarse-graining chromatin structure and quantifying effective interactions between the obtained partitions. In the previous section, we have described the interaction network of a chromosome c in terms of a MSM with transition matrix pij (Eq 2), representing the probability for a probe to reach locus j from locus i, with the steady-state distribution of probes μi in locus i (Eq 3), and with the undirected flux of probes πij between loci i and j defined in Eq 1. Considering the optimal hub set M with its associated committor probability qa (i) in the chromosome c, the effective flux of probes between any two soft-partitions a, b is given by Fab=∑i∈cqa (i) πib (11) where πib is the undirected flux of probes between the loci i and b (Eq 1) in chromosome c. In other words, Eq 11 measures the portion of flux between any locus i and hub b passing through hub a (see S2 Fig for illustration of the notion of effective interactions). Eq 11 is an exact calculation on a single chromosome. We are also interested in evaluating the effective fluxes between hub loci in different chromosomes. However, Eq 11 cannot be simply extended to the whole genome as the committor is by construction qa (i) = 0 for any locus i ∉ c. As computational limitations do not allow us to calculate the exact committor for the entire genome, a mean field formulation of Eq 11 was used to estimate the effective flux between any two partitions A, B in the genome. To this end, the effective flux between partitions, irrespective of the chromosomes to which they belong is calculated as FAB=∑i∈gθA (i) ∑j∈gπijθB (j), (12) where the summations are carried over the entire genome g, πij is the flux of probes between any pair of i and j in the genome (Eq 1 with fij the Hi-C matrix of paired-end reads of counts is now extended to the entire genome), and θA (i) the hard-partitioning committor defined in Eq 8. The rationale of Eq 12 is to efficiently, though indirectly, estimate the flux between any two partitions A, B in terms of all the intermediate pairwise fluxes πij. Within the logic of a MSM, the effective fluxes in Eq 12 serve as a measure of chromatin effective interactions. Chromosome partitions are obtained from the optimal hub set as a result of the metastability analysis upon annealing conditions. They offer a coarse-grained description of the genome as the interactions between partitions are characterized via effective interaction strengths (Eq 12). Given a genome-wide set of partitions obtained above, a putative reduced model of the major partition interactions can be constructed by directly coarse-graining the matrix of counts fij for the entire genome. The observed joint probability of interaction between two partitions A and B is P (A∩B) =∑i∈A∑j∈Bfij∑ (X, Y), X≠Y∑i∈X∑j∈Yfij, A≠B (13) where the summation in the denominator is carried out on the pairs (X, Y) of distinct partitions to ensure proper normalization. Because of the law of total probability, the probability for a partition A to be involved in any interaction other than itself is P (A) =∑Y≠AP (A∩Y) (14) which by construction adds up to one over all possible partitions A. In general, in the case of independent partitions, namely with no association between them, the relation P (A ∩ B) = P (A) P (B) would hold for the interaction probability. Therefore, to provide a measure of the degree of association between partitions, we define the following affinity as CAB=P (A∩B) P (A) P (B) (15) which is a positively defined quantity. This quantity is also known as the observed to expected ratio o/e where P (A ∩ B) and P (A) P (B) are the observed and expected probabilities, respectively. In the case of CAB > 1 where the observed probability exceeds the expected, this is interpreted as a degree of association between partitions, either a contact or functional relationship. On the contrary, if CAB ≤ 1 observed and expected probabilities either coincide or the expected probability exceeds the observed one. These situations are interpreted as either no association (between partitions CAB = 1) or dissociation (partitions repel each other for CAB < 1). Thus, high values of affinity indicate a high degree of association between partitions, suggesting the presence of active binding and/or co-localization mechanisms. Intra-chromosomal pairs show very high affinities, typically with CAB > 10, while inter-chromosomal pairs have affinities CAB < 4. In this work, we analyzed 50kbp in-situ Hi-C interaction maps obtained by Rao et al. [23] for human B lymphocyte cells (GM12878, two replicates) at both single-chromosome and whole-genome levels (GEO accession GSE63525). Three other datasets listed under the same GEO accession were also analyzed: IMR90 (lung fibroblast), HUVEC (umbilical vein endothelial cells), and HMEC (mammary epithelial cells). Epigenomic data tracks for GM12878 were obtained from the ENCODE Consortium web portal, with signal tracks for transcription factor ChIP-Seq from ENCODE/Stanford/Yale/USC/Harvard, histone ChIP-Seq from ENCODE/Broad Institute, DNase-Seq from ENCODE/OpenChrom (Duke). Z-scored fractions of epigenetic factors were calculated in order to investigate their distributions within partitions. In the single-chromosome case, for a given signal track density xf (A) of factor f in a partition A of chromosome c, the Z-scored density of factor f is: Zf (A) =xf (A) −μfσf (16) where μf and σf are the weighted mean and standard deviation of densities of factor f across partitions in the chromosome c. For the Z-score calculations on the whole-genome, the weighted mean and standard deviation across all 539 partitions were used. For the network representation of the effective interactions, the force-directed layout in Cytoscape was used [79] with the force constants parametrized as logFAB, where FAB is the effective interaction between partitions (Eqs 11 and 12). The node sizes are proportional to the partition size or Z-scored epigenetic factor density, respectively. Only partitions of size larger than 2Mbp are shown. Edge width scales with logFAB and only interactions above a certain threshold are shown. For intra-chromosomal networks, width of edges is defined according to fixed thresholds of the interaction strength at each level of hierarchy. In the whole-genome network of effective interactions, given the large number of partition pairs with a wide spread of effective interaction strengths, we classify interaction strengths into discrete levels and ignore weaker interactions. Histograms of the distribution of effective interaction strengths are plotted in S4 Fig, with intra-chromosomal (red), inter-chromosomal (green), and all (blue) interactions shown on the same axis. Layers of successively weaker interactions provide finer details to the interaction network structure (see S4 Fig): (i) Scaffold-Layer interactions are the strongest 2000 interactions, or the top 1. 35% of all interactions; (ii) Layer 1 interactions comprise the top 1. 35% to 1. 5% of all interactions, compared with the Scaffold Layer; (iii) Layer 2 interactions represent the top 1. 5% to 1. 7% of interactions; (iv) Layer 3 interactions represent the top 1. 7% to 2. 0% of interactions. In our analysis for the GM12878_primary network, we considered only the scaffold and Layer 1 interactions (the top 1. 5% of all interactions) to be significant. Before performing the MSM analysis for single chromosomes, a Gaussian Filter (GF) was employed to reduce the effects of sampling noise and systematic errors in Hi-C data: the matrices of raw interaction counts were convolved with a Gaussian kernel. With the interaction matrix at 50kbp resolution, a width parameter in the Gaussian kernel σ = 200kbp was used, truncated at 4σ. S11 Fig shows a comparison of the raw Hi-C matrices with those after the GF preprocessing on chromosomes 1,17, and 20. Performing partitioning analysis with and without GF preprocessing showed that both approaches yielded similar hub sets and partitions, and also when σ is varied within reasonable bounds, but the optimization on filtered datasets converged more rapidly. The Gaussian kernel width σ was chosen to balance between retaining structural information and computation speed: while increasing σ improved convergence rate, doing so smears out structural information in the high-resolution interaction matrices. Computation of effective interactions between partitions (Eq 12) is not affected directly by GF as the raw interaction matrices are used for obtaining the πij values. The algorithms used in this study are implemented in a freely available Python package ChromaWalker (https: //bitbucket. org/ZhenWahTan/chromawalker), built on the standard SciPy stack of libraries (NumPy, SciPy, Matplotlib, and Pandas), using a serial implementation on CPU. The run time for a full genome at 50kbp resolution, on a 3. 4GHz Intel Core i7 CPU with 8GB RAM, is approximately 1 week.
A new era in chromatin research started with the availability of Hi-C data and new experimental techniques driving improvements in data resolution enable us to achieve a deeper understanding of the chromatin structure and function, while calling, at the same time, for the development of more advanced analytical methods. Though instrumental in the analysis of Hi-C data, both model-driven polymer models and data-driven statistical approaches are always based on several assumptions and require tweaking parameters. We interpret the Hi-C frequencies of chromatin interactions in terms of pairwise contact energies, obtaining corresponding energy landscape that represents the structure and interactions in chromatin. The ruggedness of this landscape is explored by the random walk of a travelling probe, which is formalized in the framework of a Markov State Model. The multilevel energy landscape induces metastability in the Markov process, revealing the hierarchy of chromatin structural organization. Structural partitions determined by the basins in the energy landscape are, thus, naturally obtained at different levels of hierarchy without any preliminary assumptions. Effective interactions between partitions are evaluated, providing a blueprint of the whole-genome organization and functional interactions, which is further substantiated by mapping of information on gene expression regulators and different epigenetic factors. The notion of affinity between partitions complements the picture by reflecting the degrees of association between partitions, calling for the modelling of chromatin dynamics and exploring its functional modulation.
markov models
chromosome structure and function
statistics
histone modification
mathematics
epigenetics
chromatin
chromosome biology
gene expression
chromatin modification
genetic loci
probability theory
chromosomes
cell biology
genetics
biology and life sciences
physical sciences
genomics
markov processes
Exploring chromatin hierarchical organization via Markov State Modelling