Circuit topology analysis of cellular genome reveals signature motifs, conformational heterogeneity, and scaling

Subject areas: Genomics, In silico biology, Mathematical biosciences

First-order intercellular heterogeneity: cells present different percentages of topological entanglement

Here we show how CT can be applied to genomic data to reveal cell-to-cell heterogeneity in the topological content. To this end, we apply CT analysis to 3D polymer models retrieved from Hi-C data (Stevens et al., 2017). However, the dependence of the CT framework on contacts alone makes it particularly suitable and easily applicable to Hi-C maps directly. An extensive discussion of opportunities and limitations of this approach is shown in Figure S1. From a conceptual point of view, we can sub-divide CT relations into entangled and unentangled relations: the former expectedly facilitates local entanglement and globular structures, whereas the latter promotes delocalization of contacts and a linear spatial expansion. It is important to specify that here the word “entangled” is used in a broad sense, as in this paper we are not concerned with knot detection. Once the contacts in the 3D structure have been identified with a spatial cutoff (r = 1 particle, 100 kb), the topological relation between each pair of contacts can be stored into an NxN topology matrix, where N is the number of contacts. The topology matrix is fundamentally different from a Hi-C map. Hi-C maps are contact maps, showing the 3D distance between each pair of genomic loci. In the topology matrix, not the genomic loci but contact indexes are displayed along rows and columns: contact 1 (formed by contact sites i, j), contact 2 (formed by contact sites k,l), and so on. Contact sites are numbered as they appear along the chain, left end to right end. A detailed scheme on how to go from a contact map to a topology matrix is shown in Figure 1: intrachain contacts are identified and numbered (Figure 1B); from the linear diagram (Figure 1C), it is easy to recognize the three topological relations. This information can then be stored in the topology matrix (Figure 1D), which yields a reduction in dimensionality with respect to the contact map. The topology matrix of a chromosome (Figure 2A) presents a highly ordered distribution of topological relations, with entangled fractions clustering along the diagonal in characteristic elongated features, which we call L-patterns. These structures are highly conserved among all chromosomes (Figure S2) and are reminiscent of the plaid-like square patterns found in Hi-C maps (Durand et al., 2016; Lu et al., 2020). However, here we add topological information on top of geometrical information about loci proximity provided by Hi-C; the clustering of entangled relations close to the diagonal not only indicates the tendency of forming tightly packed domains with higher contact density but also informs us on the topological content within these domains. These clusters are highly entangled and assume a looping structure with a highly conserved arrangement. By a first-order analysis, we can characterize a chromosome by the ratio of its entangled fraction over the total of its contacts. Interestingly, the entangled fraction averaged over all chromosomes in a cell varies dramatically from cell to cell (Figure 2B), indicating that the entangled fraction is a good parameter to distinguish nuclei based on their topological properties. We also calculated the characteristic entangled fraction for each chromosome, averaging over all eight cells. Although variations are present (Figure 2C), it appears that intracellular heterogeneity is less dramatic than intercellular heterogeneity (Figure 2B). This suggests that a higher order analysis is needed in order to detect intracellular variations, such as shown below in the local L-loop analysis. Moreover, this indicates that the topology of single chromosomes is most likely widely affected by the overall topological status of the whole genome, which can be quantified in terms of entangled fraction.

Topological information can be coupled to structural information to gain insight into biological mechanisms

The information about entanglement can be traced back to chromosome structure, to gain local topological insight. This insight can then be coupled with biological data, in order to trace a correlation between topological arrangement and chromosomal function. A couple of easy mathematical operations are needed in order to project topological information onto chromosome sequence. As previously mentioned, the topology matrix has size equal to NxN, where N is the number of contacts. Therefore, each row in the matrix corresponds to one specific contact in the structure, and scanning that row, we can access all topological relations that the contact experiences. If we want to know the overall entanglement that the contact experiences, we can sum all entangled relations appearing on that specific row. This information can intuitively be projected onto a contact map, which represents a very common visualization technique for Hi-C data. The result of this procedure is shown in Figure 3A as a heatmap: each pair of contact sites on the two sides of the diagonal is assigned a color depending on the number of entangled relations it experiences in the topology matrix. Fluctuations in the amount of topological entanglement are now visible on the contact maps, allowing us also to compare the pattern thus created among different cells (Figure 3B). For ease of visualization, we zoom on a small portion of chromosome 19, corresponding to 10 Mb (Figure 3B). Once again, we detect a great heterogeneity in the entanglement pattern among the various cells. However, some general features can be identified in most cases, such as a gradual increase in entangled fraction along the first 10 Mb of sequence. Sometimes this increase is not monotonic, as is for example the case in cell 5 and cell 7. Processing these data, we can take a step further and reduce the dimensionality, going from a 2D contact map to a 1D trace of the entangled fraction to be mapped along the chromosome sequence. The procedure is again comparable to what was described for the topology matrix; each row (or column) in the contact map corresponds to all the connections a single contact site experiences, with the addition of a certain weight given by the entanglement value of each contact (represented by the color scheme in the heatmap). Therefore, summing all the values on one row provides a measure of the overall entanglement experienced by each contact site (100 kb of sequence, in this case). This easy matrix operation allows us to map the level of entanglement onto chromosome sequence, as exemplified by the plot on top of the contact map in Figure 3A. Despite its simplicity, this quick procedure yields important practical implications. One could consider to couple this parameter to biological information concerning genetic sequence and look for correlations between entanglement and genetic activity. As a proof of concept of the relevance of this procedure, we tried coupling the entanglement trace with gene abundance data measured by nuclear RNA-seq; the data were provided by Stevens et al. (2017) as supplementary information. In order to provide a good signal-to-noise ratio for this proof of concept, we utilized Hi-C population data; after applying a threshold for excluding contacts with low likelihood (those contacts with counts lower than 250 in the Hi-C map), the CT pipeline thus far described could easily be generalized and implemented for this kind of data. Entanglement data were binned in 1 Mb segments and compared with the RNA-seq data for the corresponding Mb. Figure 3C shows the result of this analysis for chromosomes 4 and 6; the plots for all other chromosomes are displayed in Figure S3. The Pearson correlation coefficient highlights a statistically significant positive correlation between the two parameters, indicating that high entangled fraction might correlate with a higher activity in gene expression. This trend is confirmed for all chromosomes except two, chromosome 10 and 17 (Table 1), which show no correlation. To get a cleaner reading, we set a threshold on gene abundance as well: reads that were lower than 50 per Mb were set to zero. However, overall the correlation coefficients did not change significantly; correlations calculated without threshold for gene abundance are shown in Table S1. Similarly, correlations with different thresholds for Hi-C counts are displayed in Tables S2, S3, and S4 (thresholds set respectively to 100, 200, and 300 counts). The overall trends are not dependent on threshold choices. These results suggest a biological relevance associated with entanglement, as well as confirming the versatility of the method, which can be applied to any biological data from genetic sequence.

Second-order intercellular heterogeneity: concerted relations and clustering coefficient detect loop extrusion activity

Three-contact topological arrangement appears to be a very powerful tool for genome topology classification, allowing to identify distinct topological states with different physical characteristics. A very handy and well-established mathematical tool for studying higher order topological arrangement is the network clustering coefficient (Kaiser, 2008; Watts and Strogatz, 1998). The clustering coefficient is a measure of the cliquishness of a “neighborhood” in the network: it reflects the extent to which the neighbors of a node (its connections) are also connected with each other. For example, in a social network the clustering coefficient would indicate how many of a person's friends are friends with each other. Let us consider a node n with k_n neighbors; the maximum number of edges that can exist among these neighbors is K_n = k_n(k_n−1)/2 (Watts and Strogatz, 1998). The local clustering coefficient C_n represents the fraction of these edges that actually exists: C_n = e_n/K_n, where e_n is the number of existing edges in the neighborhood of n. The clustering coefficient was already applied successfully to Hi-C interaction data in combination with deep learning for the identification of sub-compartments (Ashoor et al., 2020). Here, we exploit the concept of the average clustering coefficient, which means that we consider the average of C_n over all N nodes in the network. We show how the clustering coefficient can detect cellular heterogeneity and exploit CT to explain what is the structural significance of such a parameter.

The interaction network is created by treating contact sites as nodes and contacts as edges. Therefore, coordinates in the polymer that do not participate in any contact are discarded and not considered in the network representation. Now, if we only considered S, P, and X relations, we would obtain a disconnected graph, where edges have only pairwise connections. Such a network would yield a null average clustering coefficient. Therefore, the only topological relations that contribute to the clustering coefficient are concerted relations (Figure 1A), as each node needs to be in contact with at least two other nodes. By looking at the case studies in Figure 4, it becomes apparent that only one possible configuration contributes to the average clustering coefficient of the network, formed by two CS contacts enveloped by a CP contact, which we will call trefoil configuration. If the network were composed exclusively by trefoils in series, it would yield a clustering coefficient of 1 (Figures 4A and 4B); on the other hand, if the network were completely devoid of this configuration, the clustering coefficient would be zero (Figures 4C and 4D). The presence of trefoil configurations embedded in structures with other topological relations modulates the clustering coefficient between 0 and 1 (Figures 4E–4H).

Sanborn et al. (2015) previously detected such “cliques” (trefoils) in networks retrieved from population Hi-C data representing an interphase mammalian genome. By combining experimental and computational assays, they were able to identify cohesin-mediated loop extrusion as the mechanism giving rise to this particular network configuration. A model for loop formation based on extrusion complexes associating to DNA was suggested (Alipour and Marko, 2012; Nasmyth, 2001) and later extended with the assumption that each subunit of the extrusion complex recognizes an appropriately oriented motif on a DNA strand, such as a CTCF motif (Sanborn et al., 2015). Such motifs work as “anchors” for the loops, halting the extrusion process whenever the extrusion complex reaches, on the two ends of the loop, a convergent motif pattern (the motifs point toward each other). The trefoil diagram corresponds to two loops coming together at their ends in the 3D space (Figure 5A). Therefore, this picture would correspond to clique loci colocalizing in the same area in the nucleus, giving rise to chromatin rosettes (Figure 5B). Sanborn et al. hypothesized also the formation of a third loop between convergent CTCF motifs, arising from loci A and C. They argued, however, that most likely two consequent loops (A-B, B-C, Figure 5B) tend to form in some cells, whereas the outer loop (A-C) in others. Because we define contacts based solely on spatial cutoffs, we cannot discriminate between an A-C loop formed by loop extrusion complexes, from one created by accidental spatial proximity. However, a cell where mostly (A-B, B-C) loops tend to form would yield in our analysis very different clustering coefficient values than one where mostly A-C loops form, as the latter would not give rise to trefoil configurations.

In order to estimate the prominence of such structures in single cells, we calculated the average clustering coefficient for every chromosome and averaged it for each cell (Figure 5C); astonishingly, clustering values for each cell divide the population neatly into two subgroups, 4 cells with high and 4 with low clustering (approximatively above and below 0.02). There is room for speculation as to what this sharp division might mean. It is possible that we are indeed observing cells where mostly A-B and B-C loops arise (high clustering), as opposed to cells where A-C loops arise, confirming Sanborn et al.’s (2015) hypothesis concerning the two processes being mutually exclusive. This mutual exclusion might explain why we get such stark sub-division among clustering coefficient values among cells.

These two populations display distinct topological properties. We find a weak but statistically relevant correlation between average clustering and entangled fraction (r = 0.31, p = 0.01) for high clustering cells, whereas there is no correlation for the low clustering group (Figure 5D). Similarly, the gyration radius of chromosomes correlates negatively with entangled fraction for high clustering cells (r = −0.45, p = 0.00), while again, there is no correlation for the low clustering cells (Figure 5E). This result is somewhat surprising, as some degree of correlation between geometrical and topological properties would be expected, being the nucleus a tightly packed constrained space (Kinney et al., 2018; Satarifard et al., 2017). These findings suggest that loop extrusion and the formation of chromatin rosettes might be related to wider changes in the overall genomic physical and geometrical properties.

Higher-order intercellular heterogeneity: fractal dimension of the topology matrix quantifies the topological state of the chromosome

Here, we search for fractality in the topological domain and show that the topological fractal dimension D can be used as a parameter to quantify the topological distance between chromosomes. Moreover, D can highlight the overall topological features of genomes and how these features differ from those of randomly generated folded chains. Unlike the approaches presented in the previous sections, the fractal dimension D accounts for both concerted and nonconcerted relations, as well as more complex topological features such as L-patterns. Fractal properties of genome structure have long been suggested; the fractal globule is one of the most used physical models employed to describe genomic 3D arrangement (Iashina and Grigoriev, 2019; Mirny, 2011; Rao et al., 2014). Moreover, the fractal dimension of coding and noncoding sequences along the linear human genome sequence was calculated (Garte, 2004), revealing ranges of lengths where this is constant and remarkably consistent among chromosomes. Here, we take a leap from geometrical properties, over which the fractal dimension is typically calculated, and consider the topological space exclusively. The topology matrix is not dependent on geometrical parameters such as distance between the contacts and is therefore particularly suitable for the identification of recurring structural patterns. We can confront a typical matrix derived from a chromosome (Figure 6A), to one derived from a randomly generated chain (Figure 6B), with comparable number of contacts. The random chain was created by performing a random walk inside a box, where the number of steps (in units of 100 kb particles) corresponds to the typical length and the size of the box to the typical size of a chromosome.

The two matrices look remarkably different, especially in the number and prevalence of entangled fractions, which are confined in ordered structures along the diagonal in the case of chromosomes. This arrangement indicates a more linear development in space of genomic loops and a higher degree of segregation of contacts inside domains. In order to quantify these qualitative differences, the fractal dimension was calculated over the two-dimensional topology matrix in its Minkowski–Bouligand variant (Falconer, 2014). The calculation of the fractal dimension over subsequently bigger portions of the topology matrix revealed a constant fractal dimension (ranging from 1 to 1.7) over a wide range of length scales, revealing that the matrix elements form a self-similar fractal in this range (Figure S4). However, the fractal dimension of the chromosomes was found to correlate positively with the total number of contacts in the chain (r = 0.41, p = 7.1 × 10⁻⁸); therefore, in order to compare chromosomes with each other, we normalized the fractal dimension by the number of contacts. The dataset was divided into a high clustering and a low clustering group, as identified by the average clustering coefficient. The results are shown in Figure 6C: the two distributions for normalized fractal dimension for high and low clustering cells are not Gaussian, have equal variance, and are statistically distinct (p = 0.023). Low clustering chromosomes (mean: 113 ± 4 × 10⁻⁵) display a statistically smaller normalized fractal dimension than high clustering ones (mean: 125 ± 5 × 10⁻⁵). This result again suggests that these two groups of cells might belong to two distinct states from a topological point of view.

Even more stunning are the topological differences between chromosomes and random chains in terms of the normalized fractal dimension (Figure 6C). Ten different random extractions of 80 random chains yielded a much wider distribution of values, indicating that the topology of chromosomes occupies a very narrow range of the space of all possible configurations. The random distribution is significantly different from both chromosome distributions for each extraction (for individual p values see Tables S5 and S6). The displacement of the random chain distribution toward higher values can be explained by the tendency of random chains to form more entangled long range contacts as opposed to chromosomes. In fact, our simulated random walks in a box create equilibrium conformations of polymer globules, which are known to correspond to a collection of overlapping random walks (Emanuel et al., 2009; Lua et al., 2004). According to Flory's theorem, each chain section shows ideal chain statistics up to extensions when it reaches the surface of the globule where it is reflected. This means that for sufficiently large chemical distances between monomers (i.e. distances along the chemical backbone) their contact formation probability is random, especially independent of the chemical distance.

The average normalized fractal dimension is also, as the entangled fraction, a valuable tool to distinguish individual cells (Figure 6D), indicating that each nucleus has its own topological fingerprint. We also observe a weak negative correlation between fractal dimension and gyration radius (r = −0.29, p = 0.00, Figure 6E). Interestingly, a correlation between the average clustering and the fractal dimension can only be observed, once again, for the high clustering group (r = 0.42, p = 0.00, Figure 6F), indicating that, for this distinct topological state, there is an interplay between the abundance in trefoil configurations and the overall topological properties of the chromosome.

Cumulative CT parameter traces reveal characteristic topological length scales

The trend of entangled fraction and fractal dimension measured over subsequently larger portions of the polymer chain reveals evidence of scale invariance in the topological content. We created windows of increasingly higher number of particles, with an increment of five particles (5, 10, 15, 20, …, n-5, n) over which to calculate topological fractions and the fractal dimension. This procedure is performed starting from the left end of the chain and gradually widening the window of analysis to include the right end. Results for the same analysis carried out by processing the chain in the opposite direction can be found in supplementary information (Figure S5). As stated previously, each particle corresponds to 100 kb. This analysis highlighted some very interesting periodical features.

By gradually increasing the length scale at which we analyze the chromosome structure, we observe the periodical appearance of local maxima in the entangled fraction (Figure 7A). The tracking of local maxima allowed for the retrieval of the distribution of length scales over which these fluctuations occur (Figure 7B). These distributions are remarkably similar for all chromosomes, with peaks around 10–12 Mb (100–120 particles). There is also no significant difference in the characteristic length scale between the high and low clustering groups (Figure 7C), indicating that this might be a structural feature that is highly conserved among genomic topological states. These conclusions allow for speculations concerning the genomic hierarchical structure, namely the existence of structures with a characteristic length scale of about 10 Mb. What kind of structure in interphase chromosomes is responsible for setting such a scaling? It is possible that this characteristic length might actually still be a "shadow" of the mitotic chromosome state. It has been estimated that the relaxation time of human chromosomes is on the order of 500 years (Rosa and Everaers, 2008), far beyond the lifetime of an organism. This timeline suggests that some features of the mitotic state survive into the interphase state. It has been speculated that both the internal conformations of interphase chromosomes and the fact that each chromosome lives in its own territory can be explained by the fact that the mitotic state of the chromosomes is remembered so that they resemble a dense solution of nonconcatenated polymer rings (Rosa and Everaers, 2008; Sazer and Schiessel, 2018). It seems thus reasonable to speculate that another characteristic feature of mitotic chromosomes (Gibcus et al., 2018), namely their internal helical structure with a repeat length of about 12 Mb, by which the sister chromatids face each other, might still be remembered by the interphase chromosome in the periodic structure of its entangled fraction. This characteristic length scale might be embedded in methods for the detection of structures such as sub-compartments (Ashoor et al., 2020). Moreover, this evidence suggests that between TADs and compartments there might be a more gradual hierarchy of structures, at least from a topological point of view, rather than a sharp structural distinction (Serizay and Ahringer, 2018).

After having defined a characteristic length scale, the question arises whether there is a length scale beyond which genome topology displays scale invariance. An object displaying fractal geometry yields the same fractal dimension regardless of the scale at which the object is observed. Clearly, a chromosome is not a pure geometrical fractal, but an object in the physical world partly displaying fractal properties (Garte, 2004). Therefore, there is a finite interval where the fractal dimension of a chromosome can be considered constant (Figure 7D). We can build a trace of the fractal dimension with the same method used for the entangled fraction and define a threshold after which D can be considered constant. By tracking the threshold of a chosen chromosome (chromosome 1) in all eight cells, we find that chromosomes in high clustering cells converge sooner to fractal behavior than those in low clustering cells (Figure 7E), i.e., a smaller number of kb is needed for the chromosome to display fractal properties. These findings should not be surprising if we consider that high clustering cells have a higher contact frequency than low clustering cells (Figure 7F); the distribution of the number of contacts in these two cells is also statistically different (p = 0.0001).

Recurring motifs in chromatin looping reveal inter- and intracellular heterogeneity

Recurring patterns in the chromosome topology matrix suggest a specific multiscale structure in chromatin arrangement. We call an L-pattern a feature of the topology matrix displaying elongated parallel branches interrupted by cross stripes (Figure 8A). These patterns tend to appear as nested in the matrix, often with cross stripes appearing along the same row/column in subsequent patterns (Figure 8B). We suggest a fold model, the L-loop, which can yield such highly conserved patterns in a topology matrix (Figure 8C). The loop originates with one (or more) parallel contact, which appears as a turn in the fold (contacts 1 and 2 in Figure 8C). The number of turns defines the thickness of the corresponding L pattern in the matrix. The (one or more) turns are then followed by a smaller contact, which we call a pocket (contacts 3,6,9 in Figure 8C). This structure, turn + pocket, can be repeated a few times in a row, as we proceed along the L-loop. In the diagram, these structures appear in parallel with each other (Figure 8D) and ultimately are responsible for the number of nested L-patterns in the matrix. In this case, we have three subsequent turn + pocket motifs, yielding three nested L-patterns in the matrix in 8A. One or more cross contacts may arise, when the loop structure thus created comes into contact with other regions of the DNA strand (contact 8 in Figure 8C). These contacts form the cross stripes perpendicular to the L-patterns in the matrix. The loop develops further after the cross contact, by forming either simple pockets (as is the case of contact 9 in Figure 8C) or new sub-loops. Inspection of the matrices (Figure 8B) shows that often one L-pattern is striped by several cross contacts, indicating that the L-loop comes into contact often with other segments of DNA, effect probably due to the spatial constraints imposed by the nucleus. The total number of contacts enveloped (in parallel relation) in a turn provides the length of the L-pattern in the matrix. For example, turn 1 envelops all other contacts, resulting in a length of 8 in the outer L-pattern. To summarize, Figure 8E displays the basic build of an L-loop: a number of turns followed by a pocket and a cluster of contacts (in red in the Figure), which can contain a variable number of turn + pocket structures, cut through by possible cross contacts (represented with a dashed line). Although the topology matrix presents several small variations of the pattern suggested here (Figure S6), this representation is a useful ideal motif to interpret looping structures in the matrix.

Although these looping structures are present in all chromosomes analyzed in this study, their numbers and lengths vary substantially for both the high and low clustering group. Considering the wide range of sizes that L-loops encompass, it was beneficial to divide them into two quantiles according to their size, above and below the median (Figures 8F–8G). The high clustering group displays a long tail toward higher values, reaching a maximum of 200 contacts for the group above the median. Figure 8H shows the matrix profile, e.g. the size of L loops (number of contacts enveloped by the loop) for chromosome 1, in all eight cells. It is apparent how for the low clustering group not only the total number of L-loops is significantly lower (unsurprising, given the total number of contacts is overall less) but also the contact frequency within the loops is strikingly different.

Considering these results in light of the aforementioned link between clustering coefficient and loop extrusion activity, one might try to draw a relation between the suggested structure for the L-loop and loop extrusion mechanisms. L loops originate by one main contact (the turn), which then encompasses all other sub-loops in parallel relation. This mechanism might very well indicate the formation of contact domains by loop extrusion (Fudenberg et al., 2016), implying therefore that the original “turn” is the loop formed by the extrusion complex. Moreover, this structural insight allows us to shed a new light on the positive correlation found earlier in our analysis between local entangled fraction and gene expression. Entangled fraction along rows corresponds to the number of contacts found in an L-loop; it is known that contact propensity participates in a bidirectional correlation with gene expression/regulation (Greenwald et al., 2019). Therefore, it is not surprising that a higher number of contacts in an extruded loop might display correlative patterns with transcriptional activity.

The length of L-patterns in the topology matrix can also provide interesting local information, once the chromosome data have been sectioned according to the desired spatial resolution and with an appropriate spatial cutoff. As a proof of concept, we divided the chain into four segments of equal number of beads and computed the total number of contacts contained in L-loops for the section. An example of this analysis carried out with eight segments instead of four is shown in Figure S7. This procedure resulted in chromosome “bar codes” (Figure 9A). It is notable, from the bar code of a single chromosome belonging to different cells, how not only the overall topological properties but also the local arrangement can vary wildly. Interestingly, the bar codes for the total number of contacts clustered in L-loops do not mirror those for the total number of contacts (Figure S8), indicating that there is no direct correspondence between the number of contacts and their arrangements into L-loops. Therefore, the information is far from redundant. The information for all cells and all chromosomes can be conveniently visualized in heatmaps, to identify trends. For example, one might choose to plot where the maximum L-loop contacts appear in each chromosome (Figures 9B and 9C). At a first glance one notices a prevalence of white in the low-clustering cell map (Figure 9B) and of purple in the high clustering one (Figure 9C), indicating that a good portion of chromosomes tend to have their maximum in their last segment when in the low clustering state, whereas they reach it in the second segment when highly clustered. We can ideally divide the chromosomes based on this tendency: Chromosome 1, 2, 3, 5, 8, 10, 11, 18, and 20 reach maxima in segment 2 (purple in the heatmap) in at least one of the cells in the low clustering group, whereas all other chromosomes do not. Therefore, two distinct groups of chromosomes are created based on their local topological arrangement. Quite interestingly, these two groups display different topological properties in both clustering cell groups. Once again, we couple our analysis to a well-known concept in network topology, network connectivity (Figures 9D and 9E). The chromosome group that does not display a maximum length of L-loops in the second segment in the low clustering state has overall higher connectivity, in both cell groups. High connectivity in gene networks was previously shown to be highly informative on disease heritability (Kim et al., 2019). Therefore, local topological arrangements can be a useful predictor about the overall topological properties of the chromosome and potentially about its biological functions as well.

Limitations of the study

The future applicability of the method strongly relies on its adaptability to Hi-C contact maps, as these constitute the most common data type in the field. The pipeline described in this paper is directly applicable to single-cell contact maps almost completely, with the exception of the “three-contact arrangement” step, which involves the calculation of the average network clustering coefficient. Contact maps store information about contacts in a binary fashion (Figure S1); therefore, they might not carry enough information about contacts formed by triplets of contact sites, such as that depicted in Figure 5A. For this reason, the average clustering coefficient calculated directly over the contact maps (as opposed to the 3D model) of the chromosomes in this paper yielded 0 in all cases. This issue might be resolved by choosing a coarser binning of the contact maps but that would come with a loss in resolution. The issue, however, does not concern contact maps from bulk Hi-C, where the presence of contact counts and likelihoods makes it possible to detect three-contact arrangement and yields therefore a nonzero clustering coefficient. The extension of the pipeline to population Hi-C is possible with the substitution of the spatial cutoff with a cutoff in number of counts. However, a rigorous study on the effect of such a choice in cutoff is still missing and will be the object of future research.

The dataset used in this paper also presents some limitations. The resolution and number of cells could be improved in order to generalize the findings of the study and gain insight into the topological arrangement of sub-TAD structures.

Key resources table

Resource availability

Method details

Quantification and statistical analysis

Author contributions

AM conceived the project and supervised the research. BS carried out the study, wrote the codes, and collected the data. AG offered feedback and suggestions on mathematical procedures. AM, BS, and HS interpreted the data. BS, HS, and AM wrote the original version of the manuscript, which was finalized with inputs from all authors.

Declaration of interests

The authors declare no competing interests.

Document S1. Figures S1–S9 and Tables S1–S6

1. Alipour E., Marko J.F. Self-organization of domain structures by DNA-loop-extruding enzymes. Nucleic Acids Res. 2012;40:11202–11212. doi: 10.1093/nar/gks925. [DOI] [PMC free article] [PubMed] [Google Scholar]
2. Ashoor H., Chen X., Rosikiewicz W., Wang J., Cheng A., Wang P., Ruan Y., Li S. Graph embedding and unsupervised learning predict genomic sub-compartments from HiC chromatin interaction data. Nat. Commun. 2020;11:1173. doi: 10.1038/s41467-020-14974-x. [DOI] [PMC free article] [PubMed] [Google Scholar]
3. Brandão H.B., Paul P., van den Berg A.A., Rudner D.Z., Wang X., Mirny L.A. RNA polymerases as moving barriers to condensin loop extrusion. Proc. Natl. Acad. Sci. U S A. 2019;116:20489–20499. doi: 10.1073/pnas.1907009116. [DOI] [PMC free article] [PubMed] [Google Scholar]
4. Carrière M., Rabadán R. In: Topological Data Analysis. Abel Symposia. Baas N., Carlsson G., Quick G., Szymik M., editors. Springer; 2020. Topological data analysis of single-cell hi-C contact maps; pp. 147–162. [DOI] [Google Scholar]
5. Cavalli G., Misteli T. Functional implications of genome topology. Nat. Struct. Mol. Biol. 2013;20:290–299. doi: 10.1038/nsmb.2474. [DOI] [PMC free article] [PubMed] [Google Scholar]
6. Di Pierro M. Inner workings of gene folding. Proc. Natl. Acad. Sci. U S A. 2019;116:4774–4775. doi: 10.1073/pnas.1900875116. [DOI] [PMC free article] [PubMed] [Google Scholar]
7. Durand N.C., Robinson J.T., Shamim M.S., Machol I., Mesirov J.P., Lander E.S., Aiden E.L. Juicebox provides a visualization system for hi-C contact maps with unlimited zoom. Cell Syst. 2016;3:99–101. doi: 10.1016/j.cels.2015.07.012. [DOI] [PMC free article] [PubMed] [Google Scholar]
8. Emanuel M., Radja N.H., Henriksson A., Schiessel H. The physics behind the larger scale organization of DNA in eukaryotes. Phys. Biol. 2009;6:025008. doi: 10.1088/1478-3975/6/2/025008. [DOI] [PubMed] [Google Scholar]
9. Flavahan W.A., Drier Y., Johnstone S.E., Hemming M.L., Tarjan D.R., Hegazi E., Shareef S.J., Javed N.M., Raut C.P., Eschle B.K., et al. Altered chromosomal topology drives oncogenic programs in SDH-deficient GISTs. Nature. 2019;575:229–233. doi: 10.1038/s41586-019-1668-3. [DOI] [PMC free article] [PubMed] [Google Scholar]
10. Fudenberg G., Imakaev M., Lu C., Goloborodko A., Abdennur N., Mirny L.A. Formation of chromosomal domains by loop extrusion. Cell Rep. 2016;15:2038–2049. doi: 10.1016/j.celrep.2016.04.085. [DOI] [PMC free article] [PubMed] [Google Scholar]
11. Garte S. Fractal properties of the human genome. J. Theor. Biol. 2004;230:251–260. doi: 10.1016/j.jtbi.2004.05.015. [DOI] [PubMed] [Google Scholar]
12. Gibcus J.H., Dekker J. The hierarchy of the 3D genome. Mol. Cell. 2013;49:773–782. doi: 10.1016/j.molcel.2013.02.011. [DOI] [PMC free article] [PubMed] [Google Scholar]
13. Gibcus J.H., Samejima K., Goloborodko A., Samejima I., Naumova N., Nuebler J., Kanemaki M.T., Xie L., Paulson J.R., Earnshaw W.C., et al. A pathway for mitotic chromosome formation. Science. 2018;359:1–29. doi: 10.1126/science.aao6135. [DOI] [PMC free article] [PubMed] [Google Scholar]
14. Golovnev A., Mashaghi A. Generalized circuit topology of folded linear chains. iScience. 2020;23:101492. doi: 10.1016/j.isci.2020.101492. [DOI] [PMC free article] [PubMed] [Google Scholar]
15. Greenwald W.W., Li H., Benaglio P., Jakubosky D., Matsui H., Schmitt A., Selvaraj S., D’Antonio M., D’Antonio-Chronowska A., Smith E.N., Frazer K.A. Subtle changes in chromatin loop contact propensity are associated with differential gene regulation and expression. Nat. Commun. 2019;10:1–17. doi: 10.1038/s41467-019-08940-5. [DOI] [PMC free article] [PubMed] [Google Scholar]
16. Hadi K., Yao X., Behr J.M., Deshpande A., Xanthopoulakis C., Tian H., Kudman S., Rosiene J., Darmofal M., DeRose J., et al. Distinct classes of complex structural variation uncovered across thousands of cancer genome graphs. Cell. 2020;183:197–210.e32. doi: 10.1016/j.cell.2020.08.006. [DOI] [PMC free article] [PubMed] [Google Scholar]
17. Hagberg A.A., Schult D.A., Swart P.J. 7th Python in Science Conference (SciPy 2008) 2008. Exploring network structure, dynamics, and function using NetworkX; pp. 11–15.https://www.osti.gov/biblio/960616 [Google Scholar]
18. Heidari M., Satarifard V., Mashaghi A. Mapping a single-molecule folding process onto a topological space. Phys. Chem. Chem. Phys. 2019;21:20338–20345. doi: 10.1039/c9cp03175h. [DOI] [PubMed] [Google Scholar]
19. Heidari M., Schiessel H., Mashaghi A. Circuit topology analysis of polymer folding reactions. ACS Cent. Sci. 2020;6:839–847. doi: 10.1021/acscentsci.0c00308. [DOI] [PMC free article] [PubMed] [Google Scholar]
20. Iashina E.G., Grigoriev S.V. Large-scale structure of chromatin : a fractal globule or a logarithmic fractal? J. Exp. Theor. Phys. 2019;129:455–458. doi: 10.1134/S106377611908017X. [DOI] [Google Scholar]
21. Falconer K. Third edition. Willey; 2014. Fractal Geometry. Mathematical Foundations and Applications. [Google Scholar]
22. Kaiser M. Mean clustering coefficients: the role of isolated nodes and leafs on clustering measures for small-world networks. New J. Phys. 2008;10:083042. doi: 10.1088/1367-2630/10/8/083042. [DOI] [Google Scholar]
23. Kantidze O.L., Luzhin A.V., Nizovtseva E.V., Safina A., Valieva M.E., Golov A.K., Velichko A.K., Lyubitelev A.V., Feofanov A.V., Gurova K.V., et al. The anti-cancer drugs curaxins target spatial genome organization. Nat. Commun. 2019;10:1–11. doi: 10.1038/s41467-019-09500-7. [DOI] [PMC free article] [PubMed] [Google Scholar]
24. Kim S.S., Dai C., Hormozdiari F., van de Geijn B., Gazal S., Park Y., O’Connor L., Amariuta T., Loh P.R., Finucane H., et al. Genes with high network connectivity are enriched for disease heritability. Am. J. Hum. Genet. 2019;104:896–913. doi: 10.1016/j.ajhg.2019.03.020. [DOI] [PMC free article] [PubMed] [Google Scholar]
25. Kinney N.A., Sharakhov I.V., Onufriev A.V. Chromosome – nuclear envelope attachments affect interphase chromosome territories and entanglement. Epigenetics Epigenet. Chromatin. 2018;11:1–18. doi: 10.1186/s13072-018-0173-5. [DOI] [PMC free article] [PubMed] [Google Scholar]
26. Lengronne A., Katou Y., Yokobayashi S., Mori S., Kelly G. Cohesin relocation from sites of chromosomal loading to places of convergent transcription. Nature. 2004;430:573–578. doi: 10.1038/nature02742.Cohesin. [DOI] [PMC free article] [PubMed] [Google Scholar]
27. Lieberman-Aiden E., van Berkum N.L., Williams L., Imakaev M., Ragoczy T., Telling A., Amit I., Lajoie B.R., Sabo P.J., Dorschner M.O., et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science. 2009;326:289–293. doi: 10.1126/science.1181369. [DOI] [PMC free article] [PubMed] [Google Scholar]
28. Lu L., Liu X., Huang W.K., Giusti-Rodríguez P., Cui J., Zhang S., Xu W., Wen Z., Ma S., Rosen J.D., et al. Robust hi-C maps of enhancer-promoter interactions reveal the function of non-coding genome in neural development and diseases. Mol. Cell. 2020;79:521–534.e15. doi: 10.1016/j.molcel.2020.06.007. [DOI] [PMC free article] [PubMed] [Google Scholar]
29. Lua R., Borovinskiy A.L., Grosberg A.Y. Fractal and statistical properties of large compact polymers: a computational study. Polymer. 2004;45:717–731. doi: 10.1016/j.polymer.2003.10.073. [DOI] [Google Scholar]
30. Mashaghi A., Van Wijk R.J., Tans S.J. Circuit topology of proteins and nucleic acids. Structure. 2014;22:1227–1237. doi: 10.1016/j.str.2014.06.015. [DOI] [PubMed] [Google Scholar]
31. Mirny L.A. The fractal globule as a model of chromatin architecture in the cell. Chromosom. Res. 2011;19:37–51. doi: 10.1007/s10577-010-9177-0. [DOI] [PMC free article] [PubMed] [Google Scholar]
32. Mugler A., Tans S.J., Mashaghi A. Circuit topology of self-interacting chains: implications for folding and unfolding dynamics. Phys. Chem. Chem. Phys. 2014;16:22537–22544. doi: 10.1039/c4cp03402c. [DOI] [PubMed] [Google Scholar]
33. Nagano T., Lubling Y., Stevens T.J., Schoenfelder S., Yaffe E., Dean W., Laue E.D., Tanay A., Fraser P. Single-cell Hi-C reveals cell-to-cell variability in chromosome structure. Nature. 2013;502:59–64. doi: 10.1038/nature12593. [DOI] [PMC free article] [PubMed] [Google Scholar]
34. Nagano T., Lubling Y., Várnai C., Dudley C., Leung W., Baran Y., Mendelson Cohen N., Wingett S., Fraser P., Tanay A. Cell-cycle dynamics of chromosomal organization at single-cell resolution. Nature. 2017;547:61–67. doi: 10.1038/nature23001. [DOI] [PMC free article] [PubMed] [Google Scholar]
35. Nasmyth K. Disseminating the genome: joining, resolving, and separating sister chromatids during mitosis and meiosis. Annu. Rev. Genet. 2001;35:673–745. doi: 10.1146/annurev.genet.35.102401.091334. [DOI] [PubMed] [Google Scholar]
36. Norton H.K., Emerson D.J., Huang H., Kim J., Titus K.R., Gu S., Bassett D.S., Phillips-cremins J.E. Detecting hierarchical genome folding with network modularity. Nat. Methods. 2018;15:119–122. doi: 10.1038/nmeth.4560. [DOI] [PMC free article] [PubMed] [Google Scholar]
37. Ong C.-T., Corces V.G. CTCF: an architectural protein bridging genome topology and function. Nat. Rev. Genet. 2014;15:234–246. doi: 10.1038/nrg3663. [DOI] [PMC free article] [PubMed] [Google Scholar]
38. Orlandini E., Marenduzzo D., Michieletto D. Synergy of topoisomerase and structural-maintenance-of-chromosomes proteins creates a universal pathway to simplify genome topology. Proc. Natl. Acad. Sci. U S A. 2019;116:8149–8154. doi: 10.1073/pnas.1815394116. [DOI] [PMC free article] [PubMed] [Google Scholar]
39. Pouokam M., Cruz B., Burgess S., Segal M.R., Vazquez M., Arsuaga J. The Rabl configuration limits topological entanglement of chromosomes in budding yeast. Sci. Rep. 2019;9:1–10. doi: 10.1038/s41598-019-42967-4. [DOI] [PMC free article] [PubMed] [Google Scholar]
40. Pucéat M. Capturing chromosome conformation. Methods Mol. Biol. 2021;2157:1–7. doi: 10.1007/978-1-0716-0664-3_1. [DOI] [PubMed] [Google Scholar]
41. Rao S.S.P., Huntley M.H., Durand N.C., Stamenova E.K., Bochkov I.D., Robinson J.T., Sanborn A.L., Machol I., Omer A.D., Lander E.S., Aiden E.L. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell. 2014;159:1665–1680. doi: 10.1016/j.cell.2014.11.021. [DOI] [PMC free article] [PubMed] [Google Scholar]
42. Rosa A., Everaers R. Structure and dynamics of interphase chromosomes. PLoS Comput. Biol. 2008;4:e1000153. doi: 10.1371/journal.pcbi.1000153. [DOI] [PMC free article] [PubMed] [Google Scholar]
43. Rougier N.P. 2016. fractal-dimension.py.https://gist.github.com/rougier/e5eafc276a4e54f516ed5559df4242c0 [Google Scholar]
44. Rowley M.J., Corces V.G. Organizational principles of 3D genome architecture. Nat. Rev. Genet. 2018;19:789–800. doi: 10.1038/s41576-018-0060-8. [DOI] [PMC free article] [PubMed] [Google Scholar]
45. Sanborn A.L., Rao S.S.P., Huang S.C., Durand N.C., Huntley M.H., Jewett A.I., Bochkov I.D., Chinnappan D., Cutkosky A., Li J., et al. Chromatin extrusion explains key features of loop and domain formation in wild-type and engineered genomes. Proc. Natl. Acad. Sci. U S A. 2015;112:E6456–E6465. doi: 10.1073/pnas.1518552112. [DOI] [PMC free article] [PubMed] [Google Scholar]
46. Satarifard V., Heidari M., Mashaghi S., Tans S.J., Ejtehadi M.R., Mashaghi A. Topology of polymer chains under nanoscale confinement. Nanoscale. 2017;9:12170–12177. doi: 10.1039/C7NR04220E. [DOI] [PubMed] [Google Scholar]
47. Sazer S., Schiessel H. The biology and polymer physics underlying large-scale chromosome organization. Traffic. 2018;19:87–104. doi: 10.1111/tra.12539. [DOI] [PMC free article] [PubMed] [Google Scholar]
48. Scalvini B., Sheikhhassani V., Mashaghi A. Topological principles of protein folding. Phys. Chem. Chem. Phys. 2021;23:21316–21328. doi: 10.1039/d1cp03390e. [DOI] [PubMed] [Google Scholar]
49. Scalvini B., Sheikhhassani V., Woodard J., Aupič J., Dame R.T., Jerala R., Mashaghi A. Topology of folded molecular chains: from single biomolecules to engineered origami. Trends Chem. 2020;2:609–622. doi: 10.1016/j.trechm.2020.04.009. [DOI] [Google Scholar]
50. Schullian O., Woodard J., Tirandaz A., Mashaghi A. A circuit topology approach to categorizing changes in biomolecular structure. Front. Phys. 2020;8 doi: 10.3389/fphy.2020.00005. [DOI] [Google Scholar]
51. Serizay J., Ahringer J. ScienceDirect genome organization at different scales : nature , formation and function. Curr. Opin. Cell Biol. 2018;52:145–153. doi: 10.1016/j.ceb.2018.03.009. [DOI] [PubMed] [Google Scholar]
52. Spielmann M., Lupiáñez D.G., Mundlos S. Structural variation in the 3D genome. Nat. Rev. Genet. 2018;19:453–467. doi: 10.1038/s41576-018-0007-0. [DOI] [PubMed] [Google Scholar]
53. Stephens A.D., Banigan E.J., Marko J.F. Chromatin’s physical properties shape the nucleus and its functions. Curr. Opin. Cell Biol. 2019;58:76–84. doi: 10.1016/j.ceb.2019.02.006. [DOI] [PMC free article] [PubMed] [Google Scholar]
54. Stevens T.J., Lando D., Basu S., Atkinson L.P., Cao Y., Lee S.F., Leeb M., Wohlfahrt K.J., Boucher W., O’Shaughnessy-Kirwan A., et al. 3D structures of individual mammalian genomes studied by single-cell Hi-C. Nature. 2017;544:59–64. doi: 10.1038/nature21429. [DOI] [PMC free article] [PubMed] [Google Scholar]
55. Virtanen P., Gommers R., Oliphant T.E., Haberland M., Reddy T., Cournapeau D., Burovski E., Peterson P., Weckesser W., Bright J., et al. SciPy 1.0: fundamental algorithms for scientific computing in python. Nat. Methods. 2020;17:261–272. doi: 10.1038/s41592-019-0686-2. [DOI] [PMC free article] [PubMed] [Google Scholar]
56. Watts D.J., Strogatz S.H. Strogatz - small world network nature. Nature. 1998;393:440–442. doi: 10.1038/30918. [DOI] [PubMed] [Google Scholar]

Supplementary Materials

Data Availability Statement

• •

  This paper analyses available published data (Stevens et al., 2017). Accession numbers are listed in the key resources table.

• •

  All original codes have been stored in a Github repository and uploaded to Zenodo: https://zenodo.org/badge/latestdoi/428377249. All codes are publicly available at the date of publication. DOIs are listed in the key resources table.

• •

  Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon reasonable request.