Exploring 3D community inconsistency in human chromosome contact networks

The Hi-C method is an experimental technique that explores the 3D structure of genomes [1]. The procedure binds together spatially interacting DNA sequence pairs, which may be several megabases apart on the linear genome. If using a cell population, experimenters obtain a collection of 'binding events' representing the ensemble-averaged contact frequencies.

In practice, Hi-C starts with 'freezing' the chromosomal 3D configuration using formaldehyde. Formaldehyde cross-links spatially close DNA segments and proteins. After removing the cell membrane, a series of steps cut the DNA into numerous fragments, merge the free ends, and tag the cut spot with biotin. These steps are known as restriction digestion, proximity ligation, and biotin filling. The final stages purify DNA from proteins and RNA, remove all cross-links, and pull down biotin-marked fragments that next undergo massive sequencing. After sequencing, the fragments are associated with linear genomic locations and then processed into a contact map ('Hi-C map'), a matrix that illustrates the interaction frequency between all loci.

2.1.1. Transforming Hi-C data into a weighted network

One of the most popular ways to detect network communities [6, 7]—densely-connected substructures—is to maximize the objective function called modularity ⁵

$$\begin{equation} \mathcal{M} = \frac{1}{2m} \sum_{i \ne j} \left[ \left( A_{ij} - \gamma P_{ij} \right) \delta(g_i,g_j) \right ] \,. \end{equation} \tag{ 1 }$$

Here, A_ij denotes the adjacency matrix elements corresponding to the interaction weights between nodes i and j ($A_{ij} = 0$ indicates no edge), and P_ij is the expected edge weight based on a priori information. The most popular choice is when considering only the overall tendency of node-node interaction $P_{ij} = k_i k_j / (2m)$, where k_i represents node i's strength (the sum of its weights) and m is a normalization constant ensuring that $-1 \leqslant \mathcal{M} \leqslant 1$. Finally, g_i is the community index of node i and δ is the Kronecker delta. A key parameter in our study is the resolution parameter γ, which controls the overall community scale [13].

In principle, maximizing the modularity function with respect to all of the possible community divisions, encoded as $\{g_i \}$ in equation (1), is a mathematically well-defined deterministic concept. However, due to the computational limitation imposed by the problem, it is prohibitively difficult to find the exact solution, e.g. from the comprehensive enumeration of the network divisions. Therefore, most network community detection algorithms rely on various types of approximations or parameter restrictions. Many algorithms take a stochastic approach to sample the community partitions, just as in standard Monte Carlo [20]. One example is the Louvain-type algorithms [21, 22] we use here (detailed in section 2.3).

Although stochastic approaches like Louvain have been successful in terms of speed and accuracy in many community detection applications, their stochastic nature may produce multiple results that sometimes include inconsistent elements. Researchers tend to work around this inconsistency [23, 24] by choosing the most consistent, or reproducible, network partition. However, one of the authors of this paper has turned this inconsistency into an advantage, using it to probe network structural information [15, 16] at both global and local levels (figures 1 and 2). In particular, by studying inconsistency measures one may pinpoint scale regimes or specific node collections that are the most statistically reliable (at the global level) or flexible (at the local level). For a detailed theoretical framework, we defer to [16]. But below, we remind the reader of the essential parts used in this analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One metric we study is partition inconsistency (PaI). It quantifies the global degree of inconsistency among community partitions in the entire network (figure 1). PaI is based on a recently developed similarity measure called element-entric similarity [25], which overcomes notable biases, for example, skewed cluster sizes in normalized mutual information (and other conventional similarity measures). The similarity measure $S_{\alpha\beta}$ between community configurations α and β is defined as

$$\begin{equation} S_{\alpha\beta} = \frac{1}{N}\sum_{i = 1}^{N}\left(1 - \frac{1}{2d}\sum_{j = 1}^{N}|f_{ij}^{\,\alpha} - f_{ij}^{\,\beta}| \right) \,, \end{equation} \tag{ 2 }$$

where N is the number of nodes and $f_{ij}^{\,\alpha}$ is node j's relative importance in the stationary state of the personalized PageRank algorithm (PPR) [26] starting from the node i in configuration α. Also, d is a damping factor we set to 0.9 (default value in [25]).

Based on equation (2), we define PaI value as

$$\begin{equation} \mathrm{PaI} = \left(\sum_{\alpha = 1}^{\mathcal{C}}\sum_{\beta = 1}^{\mathcal{C}}p_\alpha p_\beta S_{\alpha\beta}\right)^{-1} \,, \end{equation} \tag{ 3 }$$

which corresponds to the average similarity between all configuration pairs. The terms form a weighted sum of $S_{\alpha\beta}$ by the proportion p_α and p_β of each configuration pair α and β, respectively, that appear in an ensemble of community detection results composed of $\mathcal{C}$ configurations. The PaI value $(\geqslant\!\!1)$ indicates the effective number of independent configurations. A small (large) PaI value represents more consistent (inconsistent) regimes, respectively. Using PaI, one extracts the most statistically reliable ranges of community scales by focusing on the (local) minima of PaI, in particular, alongside another meaningful evidence of stable communities: the number of communities stays flat at a specific integer value (illustrated in figure 1(b)).

While PaI describes the network's global inconsistency, we use another metric, membership inconsistency (MeI), to quantify local (individual-node) inconsistencies (figure 2). The MeI measure for each node i is based on the Jaccard index defined as

$$\begin{equation} J_{i\alpha \beta} = \frac{|\psi_{i\alpha} \cap \psi_{i\beta}|}{|\psi_{i\alpha} \cup \psi_{i\beta}|} \,, \end{equation} \tag{ 4 }$$

where $\psi_{i\alpha}$ indicates the set of nodes belonging to the same community as node i in configuration α. Then, the MeI measure is defined as

$$\begin{equation} \Psi_{i} = \left( \sum_{\alpha = 1}^{\mathcal{C}}\sum_{\beta = 1}^{\mathcal{C}} p(\psi_{i\alpha})p(\psi_{i\beta}) J_{i\alpha \beta} \right)^{-1} \,, \end{equation} \tag{ 5 }$$

where $p(\psi_{i\alpha})$ is the relative frequency of appearance $\psi_{i\alpha}$ in the ensemble. MeI represents the effective number of independent communities for a specific node across different community configurations. As shown in figure 2(b), the MeI values properly detect the functionally flexible or 'bridge' nodes participating in different modules ⁶ . In section 3, we use PaI and MeI to study global and local community inconsistencies of Hi-C maps and relate to these metrics to other biological data.

The stochastic community detection method we utilize throughout our work is version 2.1 of GenLouvain [22] (see https://github.com/GenLouvain/GenLouvain for the latest version). GenLouvain is a variant of the celebrated Louvain algorithm [21], which is one of the most widely used algorithms and is popular due to its speed and established packages in various programming languages.

Starting from single-node communities, the algorithm accepts or rejects trial merging processes based on the modularity (equation (1)) change in a greedy fashion. More specifically, each iteration consists of two steps: (i) looping over the nodes, and (ii) treating each community as a 'supernode' [27]. In the first process, each node is put into the neighbors community which results in the largest modularity increase $\Delta \mathcal{M}$ compared with the current partition. The second process transforms the network into a weighted network, where each community from the first step is represented as a supernode. The edges between supernodes are assigned weights indicating the number of edges joining nodes in the corresponding groups (note that self-edges for each supernode are allowed here to represent the internal edges). These steps are repeated until increasing $\mathcal{M}$ is no longer possible by moving a (super)node to a different community.

To determine the community stability across network scales, we run GenLouvain several times, at least 100, for each resolution parameter γ and then calculate the global (PaI) and local (MeI) inconsistency metrics [16].

We use GenLouvain to produce an ensemble of community partitions from Hi-C data for fixed scale parameters γ. However, some of these partitions seem correlated. To better understand these correlations, we use a graphical embedding technique to illustrate high-dimensional data. Specifically, we use t-SNE (t-distributed stochastic neighbor embedding) [28]. t-SNE is an algorithm that maps high-dimensional data into 2D, preserving both local and global structures. This contrasts with PCA, which preserves global structures. While PCA projects the data onto the dimensions with the highest variance, t-SNE focuses on maintaining relative distances between data points. Notably, unlike PCA, the axes in a t-SNE plot lack specific meanings and represent abstract dimensions used to retain the data structure in the reduced space. This property makes t-SNE ideal for visualizing complex data sets.

The t-SNE algorithm aggregates data points based on some distance metric. While there are several choices, we use the so-called correlation distance D, which is common for random vectors and defined as

$$\begin{equation} D = 1 - r(\mathbf{u}, \mathbf{v}), \end{equation} \tag{ 6 }$$

where $r(\mathbf{u},\mathbf{v})$ is the correlation between the vectors u and v, conventionally defined as

$$\begin{equation} r(\mathbf{u}, \mathbf{v}) = \frac{(\mathbf{u} - \bar{\mathbf{u}}) \cdot (\mathbf{v} - \bar{\mathbf{v}})} {{||(\mathbf{u} - \bar{\mathbf{u}})||}_2 {||(\mathbf{v} - \bar{\mathbf{v}})||}_2}, \end{equation} \tag{ 7 }$$

where $\bar{u}$ and $\bar{v}$ are the mean of the elements (so that $\bar{\mathbf{u}} = \bar{u} \times [1,1,1,\ldots]$) and $|| \cdots ||_2$ is the Euclidean norm. According to equation (6), D = 0 if they are perfectly correlated ($r = 1)$ and D ≈ 1 if they are uncorrelated (r ≈ 0).

In our analysis, we create these vectors from 100 GenLouvain runs. Each vector is a binary representation of the node-community membership for one node (each element is 1 or 0) depending on whether the node belongs to a specific community at a particular GenLouvain iteration. For our analysis, we used scikit-learn [29]. To reach the best visualization results, we tuned several parameters: perplexity: 20, early exaggeration: 8, initialization: random, and the number of iterations: 1356.

In results (section 3), we analyze the inconsistency measures PaI and MeI in terms of chromatin states derived from an established chromatin division [30] that we downloaded from ENCODE [31] (GM12878, Accession: wgEncodeEH000784). This data set constitutes a list of start and stop positions associated with chromatin states called peaks. These peaks result from integrating several biological data sets, e.g. ChIP-seq and RNA-seq, with a multivariate hidden Markov model (HMM). The authors [30] use 15 'HMM states' (S1–S15): active promoter (S1), weak promoter (S2), inactive/poised promoter (S3), strong enhancer (S4 and S5), weak/poised enhancer (S6 and S7), insulator (S8), transcriptional transition (S9), transcriptional elongation (S10), weakly transcribed (S11), Polycomb-repressed (S12), heterochromatin (S13), and repetitive/copy number variation (S14 and S15).

While identifying these states is automated (using chromHMM [32]), their names result from human interpretations, reflecting typical biological functions. Also, selecting 15 states is somewhat arbitrary, with the typical range spanning 5–20. Another limitation to consider is the resolution. It is approximately one nucleosome (${\sim}200$ base pairs) because the data stems from epigenetic histone modifications.

The start and stop regions for these 15 HMM do not match perfectly with the Hi-C bins. To classify every Hi-C bin into one of these HMM states, we calculate the folds of enrichment (FE) relative to a chromosome-wide average according to the following steps.

• 1.  
  Count the number of peaks k_X per bin, where $X = \mathrm{S1}, \dots, \mathrm{S15}$. Because some peaks span multiple bins, we only count the peak starts.
• 2.  
  Calculate the peak frequency's expected value using the hypergeometric test (chromosome-wide sampling without replacement). The expected number of X peaks per bin is calculated as $\bar k^{^{\prime}}_X = K_X \times (n/N)$, where, n is the number of peaks of any state in a bin, N is the total number of peaks per chromosome, and K_X is the total number of peaks for state X.
• 3.  
  Calculate the folds of enrichment $\mathrm{FE}_X$ for each HMM state X per bin by dividing the observed by the expected peak number, $\mathrm{FE}_X = k_X/\bar k^{^{\prime}}_X$.

We note each Hi-C bin can be enriched in several chromatin states. Based on enrichment, we divide Hi-C bins into five groups (A–D) if $\mathrm{FE}_X \gt 1$:

• (A)  
  Promoters: $X =$ S1 and S2.
• (B)  
  Enhancers: $X =$ S4, S5, S6, and S7.
• (C)  
  Transcribed regions: $X =$ S9, S10, and S11.
• (D)  
  Heterochromatin and other repressive states: $X =$ S3, S13, S14, and S15.
• (E)  
  Insulators: $X =$ S8

Apart from these five groups, we assign bins that are not enriched in any state to the category 'NA'.

In section 3, we visualize the MeI data using so-called boxen or letter-value plots. These contrasts regular box plots with outliers, where boxen plots better represent the actual distribution of the data.

In more detail, boxen plots are based on so-called letter values. These are recursively defined order statistics with specific depths d_i that splits the data set. These depths are calculated as:

$$\begin{equation} d_1 = (1+n)/2, \ \ d_i = (1 +\lfloor d_{i-1} \rfloor)/2, \end{equation} \tag{ 8 }$$

where n is the number of data points and $\lfloor \ldots \rfloor$ denotes the floor function. Here, d₁ corresponds to the median that splits the dataset into two halves. Next, d₂ splits these two halves to get the fourths. Splitting yet more times gives the eights, sixteenths, etc.

This recursive division is similar to a regular boxplot. First, we have the median in the middle, which splits the data points into two halves. Then, the remaining data spaces are halved again to get the quantiles. The boxes around the median then indicate the end of these quantiles. But boxen plots do not truncate after the quantiles. The remaining data spaces are divided yet again, adding more and more boxes to the plot that represents an increasingly smaller fraction of the dataset. The splitting stops after reaching some stopping criterion, typically the 95% confidence interval.

We illustrated the local inconsistency associated with a single Hi-C map in figure 3. This map depicts the number of contacts between all 100 kb DNA-segment pairs (KR normalized) in human chromosome 10. Along the diagonal, we highlight the GenLouvain-derived communities [22], at resolution parameter γ = 0.75. Squares sharing colors have the same community membership. These colors are better illustrated in the stripe above the map, showing how communities appear along the linear DNA sequence. We note that some scattered segments have the same color, which indicates that communities assemble DNA segments in 3D proximity, not only 1D adjacent neighbors. This contrasts the conventional notion of TADs, which comprise contiguous DNA stretches. To separate notations, we denote unbroken units of DNA stretches in a community as domains.

Zoom In Zoom Out Reset image size

It is essential to realize that the domains and communities in figures 3(a) and (b) represent a single configuration, or partition, of Hi-C network communities at one specific resolution parameter value γ. Since GenLouvain uses a stochastic maximization algorithm, we expect to find other partitions if running it several times on the same data set, some of which may differ substantially. To quantify this variability, we generated 100 independent network partitions and calculated the local inconsistency measure MeI [16] that quantifies how many different community configurations a single node effectively belongs to. We plot the MeI profile along chromosome 10 in figure 3(c). This profile shows that about half of the domains do not change community membership (the median value of MeI = 1.02), whereas the rest show significantly more variability (MeI ≈ 4). We also note that the MeI score is relatively uniform within each domain and that sharp MeI transitions occur near domain boundaries.

Based on previous work [16, 24], we anticipate that the MeI profile changes with the network scale. Therefore, we scanned through a wide range of community scales, extracted 3D communities, and calculated the MeI profile. We show the result from such a sweep in figure 4(a), where each MeI profile is associated with one γ value. We note that some DNA regions have low MeI scores (γ > 0.6), which indicates that nodes in those regions mostly appear in the same communities for most γ values. We indicate this as colored rectangles below the MeI profiles. But other DNA regions show the opposite behavior. These regions contain nodes that often do not appear in the same communities, which results in high and variable MeI values. Overall, the local node inconsistency grows as γ becomes larger.

Zoom In Zoom Out Reset image size

To appreciate the MeI variations from a biological perspective, we analyzed them relative to local chromatin states. As outlined in the Methods (section 2.5), we use five states and calculate the folds of enrichment for each node. We denote the chromatin states as promoters (A), enhancers (B), transcribed regions (C), heterochromatin and other repressive states (D), and insulators (E).

Below the MeI profile in figure 4(b), we show boxen plots for three γ values. The boxen or letter-value plot offers a more detailed visualization of the data distribution than a traditional box plot with outliers [33]. The boxes are determined by calculating a series of quantiles that successively divide the data set into more segments than the typical quartiles in regular boxplots (halves, fourths, eights, etc). This division provides a more detailed picture of the data distribution, especially revealing the behaviors in the tails. We defer to methods (section 2.6) for more details.

Each subplot illustrates the distribution of MeI scores associated with each chromatin group (A–E); 'NA' represents nodes not enriched in any chromatin type. Following the MeI medians (horizontal lines), we note that groups A–C have consistently higher values than the rest. In panel (c), we explore this observation more thoroughly and plot the median MeI for several γ values. The lines show that MeI grows with γ, where the A–C groups are more inconsistent than the chromosome-wide average (denoted 'all'). They also have higher MeI scores than the D–E groups that tend to follow the average. These two chromatin groups, A–C and D–E, represent the classical division into open (active) and closed (inactive) chromatin, or so-called eu- and heterochromatin.

We find similar patterns when calculating the MeI scores for all chromosomes (figure 5). We note that open chromatin nodes (A–C) are systematically more inconsistent than nodes categorized as closed chromatin and the chromosome-wide averages. However, there are significant differences in how the MeI scores change with γ. Some chromosomes, such as 1–12, show steady but wiggly growth. Others, e.g. 13–15, 20, and 22, exhibit a steep growth for small γ values that turns into a plateau. In some cases, such as chromosome 22, this plateau starts to decline, indicating that node-community memberships become more consistent.

Zoom In Zoom Out Reset image size

These MeI plots pose several intriguing biological questions that go beyond the scope of this paper. In supplementary material (figures S1–S22), we provide the complete MeI analysis like the one in figure 4.

The previous subsection analyzed the cross-scale local inconsistency measure (MeI). Here, we extend the analysis to all human chromosomes using the global inconsistency PaI instead of MeI (PaI yields one number per γ value instead of a chromosome-wide profile). The PaI score measures the effective number of independent network partitions (see section 2). By the mathematical construction of PaI [16], if there is no special scale of communities, the null-model behavior of PaI as a function of γ would be as follows. As γ gets larger from γ = 0, the PaI value starts to increase as the average number of communities increases enough to form a certain level of inconsistency (for γ = 0, PaI trivially vanishes because there cannot be any inconsistency for the single community composed of all of the nodes). On the other extreme case of $\gamma \to \infty$, each node tends to form its own singleton community, so again, there is no inconsistency, or PaI becomes zero. Therefore, if there is no particular characteristic scale of communities, the PaI curves against γ would be single-peaked ones without any nontrivial behavior such as local minima. In reality, there are characteristic scales of communities, where PaI reach its local minima [16], which indicate the most meaningful community scale. As our results show, the Hi-C communities also exhibit such characteristic scales. To better understand this metric, we revisit chromosome 10 before analyzing all human chromosomes.

We plot the PaI values for chromosome 10 in figure 6 as violet circles. When γ is small (${\lt}0.5$), we note that PaI has a plateau extending over several γ values. Such a plateau is ideal for stable community partitions (see figure 1). However, this case is trivial because the community comprises the entire network (PaI = 1, thus one effective community). Next, if γ increases above 0.5, PaI starts to fluctuate, which indicates that partitions become more variable. Notably, the growing trend stops at γ ≈ 0.6, and PaI decays to eventually reach a local minimum at γ ≈ 0.65. This local minimum represents relatively stable communities, hinted by the small effective number of independent partitions at that scale. As we increase γ above 0.65, the community structure becomes less and less stable, along with the rapidly growing number of communities (the green circles). However, asymptotically, the number of independent community ensembles grows slightly less than the number of communities per ensemble, indicating that they have different γ behaviors. For example, at γ = 0.9, there are 1.75 independent ensembles (effective), each of which is composed of 25 communities. As a final remark, it is essential to realize that the growing trend of PaI with γ does not necessarily imply the lack of intrinsic organizational scales. Instead, it indicates fuzzy scale transitions where we observe a short range of stable communities at the local PaI minimum.

Zoom In Zoom Out Reset image size

When examining the PaI curve for chromosome 10, we noticed one significant local minimum with a relatively stable community partition. Next, we ask if similar inconsistency patterns appear across all human chromosomes. To this end, we plotted PaI against γ for chromosomes 1–22 and X (except chromosome 10) in figure 7.

Zoom In Zoom Out Reset image size

We found several commonalities. First, similar to chromosome 10, most PaI curves have at least one local minimum and maximum. Some chromosomes even have two minima (e.g. chromosomes 1, 3, 9, 14, etc), which indicate multiple stable scales of communities. Also, when γ becomes large enough, the network enters the multi-community regime. Second, as γ grows, so does PaI. This growth indicates that community structures have become increasingly inconsistent. Although it is natural to observe higher inconsistency for larger numbers of communities as there are more possible combinations. As future work, it would be informative to check its scaling behavior across different chromosomes and classify chromosomes based on the functional shapes of PaI and the number of communities.

When analyzing the PaI and MeI scores across γ, we noted that the Hi-C networks exhibit relatively few independent partitions (e.g. PaI$_{\mathrm{chr10}} \lt 3$), and that each node belongs to just a few communities (median MeI $\lt 4$). This suggests that the community partitions are correlated. To better understand these correlations, we use a stochastic embedding technique called t-SNE that projects high-dimensional data clusters on a 2D plane (see section 2.4). In our case, the data set is the community membership per node over 100 GenLouvain runs.

We show the t-SNE analysis in figure 8 for three γ values, where each filled circle represents a network node (again, using chromosome 10). The closer two circles appear in the plot, the more correlated their node-community memberships are. As we increase γ, we note that node clusters split and that some circles become isolated. We interpret this as the ensemble of network partitions grows with γ and becomes increasingly dissimilar.

Zoom In Zoom Out Reset image size

While the clustering is identical in panels (a) and (b), we color-coded them differently to highlight specific features. In (a), the colors represent the local inconsistency score (MeI). We note that nodes with high MeI tend to separate from those with low MeI. We also see that the low MeI nodes have relatively stronger correlations, thereby forming more distinct clusters. Panel (a) also indicates that nodes with low MeI have similar node-community memberships.

In panel (b), the color-coding illustrates the chromatin type. To simplify the analysis, we consider two large chromatin groups—Euchromatin and Heterochromatin—instead of the five we used before (section 2.5). These groups reflect the traditional division into open and closed chromatin associated with active transcription and repression, respectively. In terms of our previous definitions, the two groups are:

Euchromatin: A, B, and C,

Heterochromatin: D.

Note that we disregarded group E ('Insulators') as it is associated with boundaries rather than long chromatin stretches, such as Eu- and Heterochromatin. Also, while most nodes belong to either Eu- or Heterochromatin, some are enriched in both types ⁷ . We call this group 'mixed'. Finally, there is yet another node group that does not enrich any of the two chromatin types ('NA').

In panel (b), we observe that community membership correlations are associated with chromatin type. For example, when γ = 0.75, Euchromatin and Mixed nodes separate from the large cluster and form new sub-groups. The nodes in these subgroups generally have high MeI scores indicating a larger variability in their community memberships.

Overall, panels (a) and (b) show a scale-dependent separation between communities associated with active and inactive DNA regions (red and blue nodes repel each other). This separation resembles the A/B compartmentalization but for small-scale 3D structures. Furthermore, the MeI score suggests that these structures have multiple independent ways to assemble if formed from Euchromatin nodes. This observation hints at higher structural variability of the accessible genome, which may reflect the dynamic nature of gene expression processes.