Identifying stable communities in Hi-C data using a multifractal null model

Chromosome contacts. We downloaded Hi-C data for the B-lymphoblastoid human cell line (GM12878) [21] from the GEO database (MAPQG0 dataset, 100 kb resolution) [29]. We only consider intra-chromosome contacts in our analysis. We interpret the data file as a weighted network (”Hi-C network”) with elements $A_{ij}$ in sparse form, where each node represents 100 kb of DNA, and the link weight is the measured contact count.

Before constructing the Hi-C network, we normalize the data using the standard Knight-Ruiz matrix balancing algorithm [30] so that the sum of all weights to a node equals unity. Note that the centromere (with no measured contact counts) must be removed before KR-normalization to avoid convergence issues. However, since the indices on the Hi-C matrix represent physical distances, we reintroduce the centromere before performing community detection (see Sec. II.3).

Chromatin states. To relate our results to different functional genomic regions, we downloaded a dataset from ENCODE [31] of 13 different chromatin types created from a multivariate hidden Markov model (denoted HMM states) [32, 33]. To reduce the number of states, we grouped the data into five effective categories: Promoters, Enhancers, Transcribed regions, Heterochromatin, and Insulators (see App. C for complete definitions in terms of the standard HMM states). Moreover, this data set allows us to calculate the chromatin state’s folds of enrichment (FE) for each Hi-C bin. We calculate enrichment relative to the chromosome-wide average (App. C).

To estimate community stability across replicates, we bootstrapped data that simulated noise from a typical Hi-C experiment. This noise comes from multiple sources [19], such as randomness of ligations between loci [34] or random contacts between loci due to fluctuations [6]. In Fig 2, we show histograms of the logged Hi-C contacts $\log(A_{ij})$ for four fixed distances in chromosome 10. We note that the distributions at all these distances are log-normal, albeit with their unique mean and standard deviation We use these empirical distributions as a proxy for experimental noise.

To generate bootstrapped Hi-maps with noise, we calculate the standard deviation $\sigma_{d}$ of the logged contacts for a fixed distance $d=|i-j|$ as

where $\langle\ldots\rangle_{d}$ denotes the average at the distance $d$. Next, we bootstrap new noisy contacts $A_{ij}^{\prime}$ for each bin using

where $\epsilon$ is the noise amplitude, and $N(0,(\epsilon\sigma_{d})^{2})$ denotes the normal distribution with mean $0$ and variance $(\epsilon\sigma_{d})^{2}$. As shown below, we tune the amplitude $\epsilon$ to not completely disrupt the chromosome’s community structure (e.g., the TADs), yet yield maps with realistic noise.

To find Hi-C network communities, or ”3D communities” [14], we use the Generalized Louvain method (GenLouvain) implemented in MATLAB [35]. GenLouvain searches for network partitions that maximize the modularity function $Q$, capturing local deviations from the expected background connectivity or null model. In its general form, the modularity is expressed as

where $A_{ij}$ are entries in the weighted adjacency (Hi-C) matrix, $m$ is the total weight of this matrix, $\gamma$ is a scale parameter setting the effective community size, $C_{i}$ is node $i$’s community assignment, and $P_{ij}^{(0)}$ is the null model. The last factor, $\delta_{ij}$, denotes the Kronecker-delta function and indicates that the sum calculates the partial sum of all local modularities within a community.

To calculate $Q$ and identify communities, one must specify a null model that captures the characteristic contact patterns in the data. The most common choice is random connections, known as the Newman-Girvan null model [26], but it does not account for structural aspects associated with long DNA polymer chain folded in 3D inside the cell nucleus [36, 37]. Previous work generalized the Newman-Girvin model to include a power-law with linear node separation $d$ as observed in Hi-C maps and predicted by theoretical polymer models [14]. However, the power-law exponent tends to change with distance [6, 28]. For example, Within TADs, the contacts decay on average as $d^{-0.75}$, whereas it is closer to $d^{-1.08}$ for longer distances. Therefore, if using a single decay function as the null hypothesis, one must introduce an arbitrary cutoff at which the decay exponent should change. To circumvent this problem, we propose another null model embracing the multifractal properties of Hi-C maps [24]. Admittedly, this model also has a parameter, but the fitting to Hi-C data is more systematic.

The model reads

where $H_{a}(i,j)$ represents the so-called Hierarchical domain model (HDM) from [24], calculated as

where

Here, $n=\lfloor\log_{2}(N)\rfloor$, $N$ is the size of the Hi-C matrix, $a$ and $b$ are fitting paramters ($b=0.5-a$). We use a matrix formulation of $H_{a}(i,j)$ and calculate the mean contact probability at each distance $d=|i-j|$, serving as $P_{ij}^{(0)}$ (see App. A). Finally, we fitted the parameter $a$ against our Hi-C data using least-squares and found $a=0.464(3)$ with minor variation between chromosomes. This value model accurately yields the decay exponent in Hi-C contact maps at all scales (Fig. 8).

Since the GenLouvin algorithm tries to optimize the modularity function in a Monte Carlo-like fashion, it often finds different local minima, producing slightly different partitions with similar quality [20] (also discussed in App. A). To reduce the number of samples required to get good statistics, we use the partition detected in the original data as the starting partition and then generate a collection of bootstrapped samples, see Fig. 1. This procedure reduces the number of required samples since we do not have to account for varying initial conditions.

To find realistic noise levels $\epsilon$ when bootstrapping our noisy Hi-C maps, we benchmarked against measured TAD stabilities between replicate experiments [16]. These experiments found the overlap of TAD boundaries between replicates to be 62%. We generated communities using GenLouvain, where we set the scale parameter $\gamma$ to achieve an effective community size $\hat{s}\sim 0.88$ Mb that best agrees with the TADs considered in [16] (see App. B). Then we generated sets of communities for $n_{\rm samp}$ bootstrapped samples for different values of $\epsilon$. To calculate the TAD boundary similarity between two sets of community boundaries $b_{i}$ and $b_{j}$ from two bootstrapped samples $i$ and $j$, we use the Jaccard index

which varies between 1 (identical) and 0 (dissimilar).

However, Eq. (7) has a few problems when using it to quantify boundary overlaps. For example, under extensive noise, the communities shrink until each node is its community. This situation yields $J_{\rm TAD}(b_{i},b_{j})\sim 1$, indicating a large overlap between community borders, but where the partitions no longer represent TAD-sized communities. To solve this problem, we multiply the Jaccard index in Eq. (7) by a ”penalty factor” containing the fraction of the number of communities in the bootstrapped sample compared to the original. This gives

where $b_{o}$ are the TAD boundaries calculated from the original dataset without noise.

To better understand how the penalty factor works, consider two illustrative examples. One is where the bootstrapped partitions contain the same number of communities as the original dataset but shifted along the chromosome. Here, the overlap is low, making the term $|b_{i}\cap b_{j}|/|b_{i}\cup b_{j}|$ small. Another case is when the bootstrapped samples have many small communities (in the extreme case, where every node represents a separate community). Then, the overlap between bootstrapped samples will be large, indicating a high $J_{\rm TAD}(b_{i},b_{j})$, but the number of communities compared to the original partition differs significantly. Therefore, the last term $||b_{i}\cup b_{j}|-|b_{o}||/|b_{i}\cup b_{j}|$ approaches $0$, leading to small $J_{\rm TAD}^{*}(b_{i},b_{j})$, suggesting poor overlap.

Finally, to provide a global measure of TAD stability in an ensemble of bootstrapped samples, we calculate the mean $J_{\rm TAD}^{*}(b_{i},b_{j})$ across all boundaries

where $n_{\rm samp}$ is the number of bootstrapped samples.

We show $J_{\rm TAD}^{*}$ over a range of noise amplitudes $\epsilon$ across all chromosomes in Fig. 3. As expected, we find perfect overlap ($J_{\rm TAD}^{*}\rightarrow 1$) when the noise level becomes low ($\epsilon\rightarrow 0$). Also, as $\epsilon$ increases, the overlap reduces below the experimental benchmark (0.62, dashed horizontal line) until reaching a plateau where each node represents separate communities. Using the 0.62-line, we may find numerically the critical noise amplitude $\epsilon^{*}$ for each chromosome. We note that $\epsilon^{*}$ varies between $\sim 0.05-0.15$. This parameter presents an analog to stability since a higher $\epsilon^{*}$ also indicates stronger resilience to noise, as the original TADs are still extractable from the dataset. From Fig. 3, we choose $\epsilon^{*}=0.11$ as the critical noise amplitude for all chromosomes (vertical dashed line).

Once we tuned the noise amplitude, we generated several bootstrapped Hi-C maps (see Sec. II.2) and extracted the parts of the community partition that stayed consistent under noise relative to the original dataset. The specific procedure reads as follows (inspired by [23]):

• •

  Get the original community partition $P_{o}$.

• •

  Generate a $n_{\rm samp}$ bootstrapped samples from the original dataset.

• •

  Get communities from each bootstrapped sample $P_{n}$, based on the starting partition $P_{0}$.

• •

  Identify the communities from the original partition $P_{o}$ that are also in the bootstrapped samples. To do this, we calculate the Jaccard indices between all community pairs in each sample and pick the ones with the largest Jaccard index.

• •

  Filter all pairs with a Jaccard index smaller than some cutoff $J_{\rm cutoff}$ and calculate the fraction of bootstrapped samples that each community in $P_{o}$ appears in.

Based on these steps,we calculate the stability fraction $S_{C_{i}}$, for community $C_{i}$ in the original partition $P_{o}$ as

Here, $\theta(\cdot)$ denotes the Heaviside step function, and $J(C_{i},C_{j})$ is the Jaccard index between two communities as formulated in Eq. (7) but not specific to TADs.

The cutoff $J_{\rm cutoff}$ is a free parameter. We set it to 1.0, meaning that we define a community $C_{i}$ as stable if all of its nodes remain in the original dataset and the bootstrapped samples. After setting this cutoff, we generated several bootstrapped Hi-C maps and calculated the stability fraction $S_{C_{i}}$ (Eq. (10)) for each community $C_{i}$ associated with the original dataset ($P_{o}$). As before, we set $\gamma$ to get the effective community size $\hat{s}\sim 0.88Mb$ for each chromosome.

We show the distributions of the stabilities $S_{C_{i}}$ as boxplots in Fig 4, one box per chromosome, sorted by the $S_{C_{i}}$ medians. By reading the chromosome numbers and the x-axis, we note no strong correlation between chromosome size and stability, as both large and small chromosomes appear on either end. To compare this result to experiments, we note that the same work that calculated the median TAD boundary similarity to be 62% also used an advanced similarity metric TADsim to measure the structural similarity of TAD sets between replicates [16]. They found this overlap to be lower, about 50% between replicates. This is close to our median stability between all chromosomes, which turns out to be $\sim 0.4$.

Until this point, we studied each chromosome separately. Below, we aggregate all the data and consider genome-wide averages. We also used $n_{\rm samp}=1000$ bootstrapped samples per chromosome to calculate the stability, which we bin to a resolution of 1%.

The boxplots and the distribution in Fig. 4 foreshadow that some of the communities have a more resilient community structure than others under noise. Here, we ask how much this stability depends on the local structure of the underlying network. In particular, we study the communities’ internal node strength. This metric differs from general node strength—the sum of all edge weights of a node—by restricting the sum to nodes within a community. Node strength is a good predictor of, for example, node centrality in a fully connected network [38] or search times to a random node [39, 40, 41, 42].

However, in our case, the stability of a community is not related to how well the nodes are connected to other parts of the network. Instead, we quantify the internal community connectedness. We define this as the sum of all edge weights within a community (including self-loops) normalized by the total node strength

where ”int” stands for internal node strength. Plotting $\hat{s}_{C_{i}}^{\rm int}$ for all communities across all chromosomes against the community stability $S_{C_{i}}$ (Eq. (10)), we find that high internal node strength suggests high community stability under noise (Fig. 5). This means that strong inter-connected communities are less affected by random noise during bootstrapping.

The previous section studied how stability is connected to network structure. Next, we perform a similar analysis using biologically relevant markers: chromatin states. In particular, we ask if our stability metric calculated on TAD-sized communities is associated with the enrichment or depletion of active or inactive genomic regions, which could suggest underlying mechanisms that form strongly connected communities. For example, there are several indications of insulation enrichment at TAD borders and facilitated promoter-enhancer interactions within TADs [43, 44].

To study this, we used data that classify genomic regions into 15 different functional states (see Tab. 2), which we grouped into five categories (defined in App. C). We plotted the folds of enrichment (FE) of these five groups against community stability for TAD-sized communities ($\hat{s}\sim 0.38$ Mb) in Fig. 6. For three of the groups associated with active genomic regions (promoters, enhancers, and transcribed regions), there is a clear positive relation between stability and FE. However, repressed regions (heterochromatin & repressive states and insulators) show no such correlation. This indicates that the active genomic regions are associated with strong stability towards noise or communities having strong internal connectivity.

Thus far, we have focused on the stability of communities of a particular size. However, by tuning the GenLouvain parameter $\gamma$ we can get communities of varying sizes. This allows us to study how communities at one scale nest with communities at larger scales. Such studies complement previous work finding that active chromatin congregates in communities with a significant hierarchical nestedness across community sizes $\hat{s}$ [15], in particular, hierarchical TADs [45, 46]. Inspired by this, we ask: do nodes in stable communities of one size stay within stable communities of a larger size? Or is community nesting related to stability?

To this end, we first calculated the community stability $S_{C}$ for different effective community sizes. We choose four commonly studied chromosome scales: TADs ($\hat{s}\sim 0.38$MB), A/B segments ($\hat{s}\sim 1.04$MB), A_1,2, B_1,2,3 structures ($\hat{s}\sim 49$MB) and A/B compartments ($\hat{s}\sim 60$MB) (see App. B). Next, we assign a node stability $s_{i}$ to each node $i$, which we choose as the community stability $S_{C}$. Thus, as $\gamma$ varies, we obtain $s_{i}$ values for each size $\hat{s}$. The change in $s_{i}$ reflects how stability changes over different network scales.

We show the change in $s_{i}$ over three different scale transitions as a boxplot in Fig. 7 (a) (green boxes). For relatively small communities (around the size of TADs), we note that the stability difference is smaller in smaller communities. To make a fair comparison, we calculated the same difference between two randomly chosen nodes (orange boxes). These boxes show that the stability difference is smaller in the actual data (green) than in the random case, at least for smaller communities. This indicates that stable communities at one size are generally supported by stable communities at a larger size, indicating the cross-scale nestedness of stable communities.

However, this relationship holds best as communities become smaller. At large sizes, the communities occupy a considerable part of the chromosome, implying that randomization does not change much since the probability of selecting two random nodes from the same community is high. This does not necessarily imply that the stability is not nested at larger sizes but rather that a much larger partition sample is required to get good statistics.

To better illustrate how nested stability changes between two select chromosome scales ($\hat{s}=0.38\mathrm{Mb}\to\hat{s}=1.04\rm{Mb}$), we plotted $s_{i}$ in an alluvial diagram (Fig. 7 (b)). To increase readability, we rounded the $s_{i}$ values to a precision of 0.2. We note that most flows from $\hat{s}=0.38\mathrm{Mb}$ to $\hat{s}=1.04\mathrm{Mb}$ go towards between nodes with similar stabilities or just above or below. This gives additional support to our conclusion that stability is nested across network scales. We also depict the chromatin state enrichment FE as bars to the left and right of the alluvial diagram (same data as in Fig. 6). We note that the FE levels for active chromatin (three leftmost bars) remain conserved across both scales for all six stability groups. The inactive chromatin shows similar conservation, but this like stem from them being uncorrelated to node stability (Fig. 6).