Systematic identification of interchromosomal interaction networks supports the existence of specialized RNA factories

Trans-C randomly walks through the Hi-C graph

Our goal is to algorithmically identify, from a given set of Hi-C data, a collection of genomic loci that exhibit strong trans interactions with each other (i.e., a “clique”). We represent the Hi-C data as a matrix, referred to as a “contact map,” in which each axis corresponds to the complete genome, and entries in the matrix represent Hi-C contact counts (Fig. 1B). In practice, the genomic axes are discretized using fixed-width bins. The bin size is thus inversely proportional to the effective resolution of the contact map. The contact map can be thought of as the adjacency matrix of a corresponding Hi-C graph, in which nodes are genomic loci (bins) and edges are weighted by the corresponding contact counts (Fig. 1C). Our goal is thus to find dense subgraphs in this Hi-C graph.

The problem of dense subgraph discovery arises in many application domains and consequently has been very widely studied. Depending on the exact formulation and the notion of density, theoretical computer science has shown that the problem complexity ranges from easily solved in linear time via a max flow algorithm (Khuller and Saha 2009) to computationally intractable (NP-hard) (Charikar 2000). Common techniques to approximate the latter case, to which our specific problem belongs, are greedy approaches, which iteratively select the best option available at the moment without guaranteeing that this strategy will bring the global optimal result (Charikar 2000), and semirandom models, which account for model errors by incorporating both adversarial and random choices in instance generation (Bhaskara et al. 2010). Trans-C approaches the discovery task of finding strongly interacting loci in trans using a random walk with restart algorithm (Fig. 1C). This general approach has been applied successfully in domains as diverse as web searching (Gibson et al. 2005), protein remote homology detection (Weston et al. 2004), and gene functional prediction (Mostafavi et al. 2008). Prior to the random walk operation, trans-C performs three preprocessing steps on the provided Hi-C contact map. First, to control for sequencing and accessibility biases, Hi-C counts are ICE-normalized (Imakaev et al. 2012). Second, the resulting matrix is processed using a binomial model to estimate interaction P-values based on an empirical null model that accounts for potential biases arising from chromosomal territorialization (i.e., small, gene-rich chromosomes generally occupying the nuclear interior and interacting more with each other than with large, gene-poor chromosomes, and vice versa) (Lieberman-Aiden et al. 2009; Bertero et al. 2019). This step allows the algorithm to focus on interactions that stand out from the noise. Third, the negative log P-values are used as weights for the network edges and subsequently refined using a “donut filter” (Rao et al. 2014) to highlight points that are local maxima. The postprocessed Hi-C interaction matrix is finally represented as a weighted graph, in which each node corresponds to a bin and the weight on each edge corresponds to the negative log P-value computed in the previous step. Trans-C then carries out a random walk with restart algorithm, which exploits global patterns of connectivity on the graph. Each walk is initiated from a randomly selected “seed” locus and moves from a node to a neighboring one probabilistically based on the weight of the edge. A parameter α controls the probability that the walk will restart at a new, randomly selected seed locus. Mathematically, as an infinite number of walks are performed, the frequency with which each node is visited converges to a stationary distribution. This can be computed analytically using the Perron–Frobenius theorem. We use the stationary distribution to obtain a ranked list of trans-interacting bins (Fig. 1D), because the most frequently visited nodes are the ones that interact most strongly with the seed loci. Highly ranked genes are most likely to be functionally related with the seed loci, and therefore, a putative clique is obtained by extracting the top ranked loci.

Trans-C uncovers the clustering of var genes critical for P. falciparum immune evasion

Having developed trans-C, we set out to test its ability to uncover known sets of loci that interact together in trans in three different organisms. First, we focused on the protozoan Plasmodium falciparum, the parasite responsible for the most lethal form of malaria. The three-dimensional organization of the P. falciparum genome is strongly associated with gene expression (Ay et al. 2014), particularly for genes involved in pathogenesis, immune evasion, and master regulation of gene expression (Bunnik et al. 2018). Among these are the variant antigenic repertoire (var) genes, a family of 60 virulence genes responsible for the antigenic variation of the parasite and evasion of the host immune system. Only a single var gene is active at a given time, the other var genes being maintained in a perinuclear cluster of heterochromatic telomeres (Fig. 2A; Duffy et al. 2017). This cluster is an excellent test case to validate the ability of trans-C to uncover a group of biologically important genes that colocalize in 3D from Hi-C data.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Trans-C identifies the var genes cluster in Plasmodium falciparum. (A) Schematic of P. falciparum in the trophozoite stage of its red blood cell life cycle, with a zoomed-in view of the nucleus highlighting its Rabl-like structure and the clustering of the var genes in a repressive heterochromatic cluster. (B) Contact heat map comparing trans contact counts among all 60 var genes versus 60 randomly selected 10 kb bins. Cis contacts are grayed out. (C) Performance evaluation of trans-C-mediated identification of var gene clustering. We plot the receiver operating characteristic (ROC) curve for the trophozoite life stage of P. falciparum. The var genes are uncovered by the random walk algorithm of trans-C with high area under the ROC curve (AUROC; 0.94). The cumulative distribution is statistically significant (P = 3 × 10⁻¹⁶⁵) from a null model of 1000 random walks performed from seeds selected randomly but with an equal or greater collective interaction strength (matched random; the line reports the average and the shaded area the 95% confidence interval). We also report the performance of a simpler greedy heuristic. (D) Visualization of the var gene–associated trans clique identified by trans-C in P. falciparum trophozoite. Nodes are color-coded by chromosome and sequentially numbered based on their relative position along each chromosome (expressed in megabases). Edges are color-coded based on trans interaction significance (cis contacts are not plotted). The seed loci for the random walk are indicated by solid black lines around the nodes.

To this end, we examined Hi-C for two stages of the P. falciparum life cycle, trophozoite and schizont, both of which are characterized by trans contacts between var genes (Ay et al. 2014). To visually highlight the var cluster (Gardner et al. 2002), we extracted the bins containing var genes and also drew 60 bins at random from the full set of genomic loci. The submatrix of trans contacts formed by the concatenation of the two sets of bins showed a striking contrast between the var and non-var loci, as anticipated (Fig. 2B). Next, we selected three var genes at random to act as seed loci and examined whether trans-C (with α = 0.5) could automatically identify the remaining 57 var gene loci. For comparison, we used a method based on a greedy heuristic that iteratively selects the bins that interact most strongly with the selected loci (Supplemental Methods). For each approach, we plotted a receiver operating characteristic (ROC) curve, in which each element is a genomic bin, labeled as positive (var gene) or negative (other loci) (Fig. 2C; Supplemental Fig. S1A). In both Plasmodium life stages, trans-C quickly found the majority of the var genes by ranking their corresponding bins highly: Of the top 50 bins, 28 contained a var gene, and all 60 var genes were recovered within the top 280 bins. Trans-C outperformed the greedy heuristic baseline, with an area under the ROC curve (AUROC) of 0.94 compared with 0.88 for the trophozoite analyses (Fig. 2C), with similar findings in schizont (Supplemental Fig. S1A). This demonstrates that a random walk approach is more suited to the task of identifying trans cliques even in the context of a clear example.

As a negative control, we run trans-C starting from three seed loci that were randomly reselected until a trio could be identified so that its trans subnetwork showed the same or greater total interaction strength as the one for the three var genes previously used as seeds. We repeated this procedure, which henceforth will be referred to as “matched random control,” for a total of 1000 times in order to estimate empirical P-values for the var gene–associated subnetwork: This proved extremely significant (P = 3 × 10⁻¹⁶⁵ and P = 3 × 10⁻¹⁷¹ for trophozoite and schizont, respectively) (Fig. 2C; Supplemental Fig. S1A). Visualization of the var gene–associated trans interaction subnetworks for the top 40 loci ranked by trans-C in both plasmodium life cycle stages showcased the intricacy of highly significant contacts (Fig. 2D; Supplemental Fig. S1B).

Identification of Greek islands regulating the expression of mouse olfactory receptor genes

To further validate trans-C, we turned to the mouse and its larger diploid nuclear genome. In mouse olfactory sensory neurons (mOSNs), chromatin regions associated to olfactory receptor gene clusters from 18 chromosomes make specific and robust interchromosomal contacts that increase in strength with differentiation (Lomvardas et al. 2006; Markenscoff-Papadimitriou et al. 2014; Monahan et al. 2017). These contacts are orchestrated by intergenic olfactory receptor enhancers that form a multichromosomal superenhancer driving the monoallelic and stochastic expression of a single mouse olfactory receptor gene (Fig. 3A; Monahan et al. 2019). The mOSN-specific trans contacts are arguably the strongest trans contacts in a mammalian genome known to date. The regions involved in such interactions were dubbed “Greek islands,” because they are sprinkled across the chromosomes as the tiny islands are in the Mediterranean Sea. Importantly, in horizontal basal cells (HBCs), the quiescent stem cell progenitors of mOSNs, these interchromosomal contacts are absent.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Trans-C identifies the Greek island cluster in mouse olfactory sensory neurons (mOSNs). (A) Schematic of trans contacts in a mOSN. The Greek islands form a multienhancer hub that is segregated from the inactive olfactory receptor (OR) genes. (B) Performance evaluation of trans-C-mediated identification of Greek island clustering. We plot the ROC curve for α = 0.5 in mOSNs versus their progenitors (horizontal basal cells [HBCs]). Trans-C correctly identifies Greek island clustering specifically in mOSNs in a way that is statistically significant (P = 6 × 10⁻¹⁹⁴) versus a matched random seed null model (average and 95% confidence interval from 1000 runs). (C) Aggregated heatmap of trans contacts among the top 60 loci selected by trans-C in mOSNs. Each square in the grid represents an average 250 kb bin in a Hi-C matrix of 21 × 21 bins centered at each interacting pair of loci (reference). The exhibited spot-like structure highlights the highly specific nature of the interchromosomal interactions of the Greek islands. (D) Visualization of the Greek island-associated trans clique identified in mOSNs by trans-C, showcasing the increased significance of loci interactions after differentiation of HBCs, plotted as described for Figure 2D.

We applied trans-C to mOSN Hi-C data (Monahan et al. 2019), randomly selecting five Greek islands from the previously reported list of 63 to use as seeds in order to measure the ability of trans-C to uncover the remaining 58. Besides running a matched random control, as a biological negative control we used HBC Hi-C data. As expected, trans-C successfully found the Greek islands in mOSNs (AUROC = 0.93; P = 6 × 10⁻¹⁹⁴) (Fig. 3B; Supplemental Table S1), although it failed to do so effectively in HBCs (AUROC = 0.71) (Supplemental Fig. S2A). At a false-positive rate of 10%, 95% of known Greek islands were identified in mOSNs, although we speculate that some of the false positives may actually represent previously unidentified Greek islands.

To visually verify whether trans-C detected specific interchromosomal contacts, we selected the top 60 predicted bins from the ranked stationary distribution (30% of which are known Greek islands). For each pair of loci from this set of 60, we extracted from the Hi-C data a 21-by-21 matrix centered at their interaction and then averaged these matrices (Fig. 3C). The resulting contact heatmap exhibited very strong punctuated signal in the middle, suggesting that the top 60 loci ranked by trans-C form specific interactions that are not driven by larger, nonspecific “neighborhood” features. Visualization of the Greek island–associated trans interaction subnetwork for the top 40 loci ranked by trans-C in mOSNs revealed a very dense network that greatly increases in significance when HBCs differentiate in mOSNs (Fig. 3D). Collectively, trans-C efficiently pinpoints trans cliques even in a complex eukaryotic genome.

We also used the Greek island data set in mOSNs to evaluate how strongly the behavior of trans-C depends on its primary parameter, the random walk restart probability α. We varied α between zero no restart to one (restart after every step) in small increments. We observed that the performance of trans-C was stable in the range [0.3, 0.7], whereas it deteriorated significantly in the two extremes when it approached zero or one (Supplemental Fig. S2B). This behavior is expected theoretically: When α is close to zero, the random walk restarts infrequently, and so, its stationary distribution becomes less dependent on the seeding bins and is mostly determined by the topology of the Hi-C graph. At the extreme, when α = 0 the walk is “memoryless” and entirely independent of the starting seed loci. On the other hand, when α is close to one there is little or no exploration along the graph. In this setting, the Hi-C data are essentially ignored, and consequently, no discoveries can be made. As an additional benchmarking, we evaluated the performance of trans-C with respect to the number of trans reads in the input Hi-C matrix. For this, we run trans-C with α = 0.5 using 100% of the trans contacts (about 436 million) versus decreasing subsamples down to about 87 million (Supplemental Fig. S2C). Trans-C maintained a comparable performance when 80% of the matrix was used, and an acceptable performance at 60% subsampling, whereas 40% and 20% of contacts proved insufficient, as could be anticipated. Lastly, we rebenchmarked trans-C against the greedy heuristic: Also in the context of Greek island discovery in mOSNs, our algorithm delivered a larger AUROC (0.93 vs. 0.92) (Supplemental Fig. S2D).

Dissecting the RBM20 splicing factory during cardiomyocyte differentiation

We next sought to explore the sensitivity of trans-C in a more challenging model in the human genome. We previously identified a network of gene loci that increase their association interchromosomally during cardiac development of human pluripotent stem cells (hPSCs) and are targets of the muscle-specific splicing factor RBM20 (Fig. 4A). Functional experiments indicated that the main RBM20 target, the large TTN pre-mRNA (which contains over 100 binding sites for RBM20), nucleates RBM20 foci. Secondary RBM20 targets interact in trans with TTN at RBM20 foci, which maximizes the efficiency of their alternative splicing (Bertero et al. 2019). We therefore dubbed the network a “trans-interacting chromatin domain” (TID) and the resulting structure a “splicing factory.” Of note, however, the cumulative interaction score of the TID calculated from shallow Hi-C data (about 90 million contacts) was only modestly enriched compared with a null model (P = 0.05). Thus, these interactions are less easily detected by Hi-C and are likely to be much more transient in nature compared with those involving the Greek islands.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Trans-C identifies the RBM20 splicing factory in human cardiomyocytes (CMs). (A) Schematic of the RBM20 splicing factory, a muscle-specific interchromosomal structure organized by the TTN pre-mRNA. This pre-mRNA binds to more than 100 copies of RBM20 and nucleates foci that engage with other RBM20 targets to promote their alternative splicing (blue arrows). (B) Performance evaluation of trans-C in uncovering the RBM20 splicing factory in early (day 15) versus late (day 80) CMs differentiated from human pluripotent stem cells (hPSCs; also analyzed as “day 0” baseline control). Results for late CMs are statistically significant (P = 5 × 10⁻¹²²) versus a matched random seed null model (average and 95% confidence interval from 1000 runs). Seed loci and ROC curves are based on a list of established RBM20 targets (Bertero et al. 2019). (C) Similar to B, but seed loci and ROC curves are based on loci directly bound by RBM20 as determined by eCLIP; P = 4 × 10⁻¹²⁰. (D) Aggregated heatmap of trans contacts between the top 60 loci selected by trans-C in late CMs starting from eCLIP data. Each square in the grid represents an average 100 kb bin in a Hi-C matrix of 41 × 41 bins centered at each interacting pair of loci extracted from the Hi-C data (reference). The denser region in the middle reveals the specific nature of the trans interactions at the RBM20 splicing factory. (E) Visualization of the RBM20-associated trans clique identified by trans-C in late CMs starting from eCLIP data, showcasing the increased significance of loci interactions during hPSC differentiation and CM maturation, plotted as described for Figure 2D.

We set out to test whether trans-C would reidentify the RBM20 TID in an independent, more deeply sequenced Hi-C data set of hPSC differentiation into the cardiac lineage (about 3 billion read pairs per sample) (Zhang et al. 2019). Besides various progenitors and early hPSC-derived cardiomyocytes (hPSC-CMs), this data set also contains late hPSC-CMs obtained after 80 days of in vitro differentiation. Moreover, older hPSC-CMs were FACS-purified using an expression reporter for the mature marker ventricular myosin light chain 2 (MLC-2v; MYL2 gene). We first attempted to recover the trans network of 16 RBM20 target genes from our original report, using five of them (TTN, CACNA1C, CAMK2D, KCNIP2, CAMK2G) as seeds for trans-C. Figure 4B shows the ROC curve for day 0 (hPSCs), 15 (early CMs), and 80 (late CMs). The best performance was achieved using Hi-C data from day 80 (AUROC = 0.84, P = 5 × 10⁻¹²²) (Supplemental Table S2); second was day 15 (AUROC = 0.78, P = 2 × 10⁻¹⁰⁵) (Supplemental Fig. S3B); and last was day 0 (AUROC = 0.75, P = 6 × 10⁻¹⁰¹) (Supplemental Fig. S3A). The improvement in ROC area as differentiation advances is in line with the important role of RBM20 in cardiac maturation (Guo et al. 2012). RBM20 is not expressed at day 0, moderately expressed at day 15, and maximally expressed at day 80. We note, however, that the performance at day 0 was better than random, suggesting that some structure that brings the loci close together is present even in undifferentiated cells. Benchmarking of trans-C against a greedy heuristic demonstrated a strong increase in performance for late CMs (AUROC 0.84 vs. 0.72) (Supplemental Fig. S3C), highlighting the advantages of the approach particularly to find trans cliques that do not stand out strongly from the noise.

Encouraged by these results, we decided to use trans-C to expand our knowledge of the RBM20 TID. Our original list of 16 genes was not the result of an unbiased search but rather reflected our prior knowledge of RBM20 biology: These 16 genes were the known splicing targets of RBM20 in both human and rat hearts that were also upregulated in hPSC-CMs. As an alternative strategy to identify genes involved in the RBM20 TID in an unbiased fashion, we hypothesized that such genes would encode for transcripts most strongly bound by RBM20 and thus enriched within the splicing factory. To test this hypothesis, we leveraged a recent data set that measured RBM20 binding to RNAs using enhanced UV cross-linking and immunoprecipitation (eCLIP) (Van Nostrand et al. 2016).

We used RBM20 eCLIP data from hPSC-CMs (Fenix et al. 2021) and counted the number of peaks that fall in each genomic bin (mapping RNAs to the encoding DNA loci). We selected the five bins with the most eCLIP peaks, which contained the genes TTN, SLC8A1, OBSCN, NEAT1, and LBD3. Using these as seed loci, we ran trans-C with α = 0.5 on Hi-C matrices from differentiating hPSC-CMs (Zhang et al. 2019). Our goal was to test whether trans-C would uncover the remaining 202 bins with eCLIP peaks. We note that this experimental setup is very different from the previous one. Here, we used Hi-C data to find binding sites in an orthogonal eCLIP data set. Moreover, only a single seed locus, TTN, was shared between this analysis and the one reported in Figure 4B. The resulting ROC curves (Fig. 4C) show the same trend: the best performance was at day 80 (AUROC 0.82; P = 4 × 10⁻¹²⁰) (Supplemental Table S3), the second at day 15 (AUROC 0.79; P = 1 × 10⁻¹⁰⁶) (Supplemental Fig. S3E), and the last at day 0 (AUROC 0.72; P = 4 × 10⁻⁹⁸) (Supplemental Fig. S3D), consistent with biological expectations.

Next, we performed a second performance recall analysis in the same trans subnetwork identified by trans-C from RBM20 eCLIP data, but in which we restricted the list of RBM20 targets to those whose RNA is bound by RBM20 on at least three sites and is differentially spliced in hPSC-CMs with RBM20 knocked out (Fenix et al. 2021). The resulting list of 45 high-confidence RBM20 targets was efficiently recovered in day 80 hPSC-CMs, with AUROC = 0.84 and a P-value of 2 × 10⁻¹²⁵ (Supplemental Fig. S3F), a performance improvement compared with the full list of RBM20 bound loci.

Similarly to our observation for trans contacts between the Greek islands (Fig. 3C), the aggregated contact frequency heatmap for loci involved in the RBM20-associated trans interaction subnetworks identified from RBM20 eCLIP data showed a clear punctuated pattern, supporting the spatial specificity of these interactions (Fig. 4D). Visualization of this subnetwork showed that it is quite dense and that it clearly increases in significance when hPSC differentiate in CMs, and even more when CMs mature (Fig. 4E). Similar results were obtained for the subnetwork identified from established RBM20 targets (Supplemental Fig. S3G).

Because the ENCODE Project generated a large number of Hi-C matrices for the left ventricle (LV) (The ENCODE Project Consortium 2012), we used this model to both evaluate the reproducibility of trans-C and determine whether the RBM20 splicing factory could be identified in adult, fully mature cardiomyocytes. We ran trans-C starting from the same five seed bins prioritized using RBM20 eCLIP data, and compared the resulting list of ranked bins: The Pearson's correlation for analyses in the 10 biologically independent LV samples was very high (range 0.77–0.85) (Supplemental Fig. S3H), whereas negative control analyses in noncardiac samples showed a low correlation (range 0.50–0.70) (Supplemental Fig. S3H). The AUROC for recovering high-confidence RBM20-bound mRNAs in LV Hi-C data was high (range 0.73–0.79), further supporting the existence of a measurable RBM20-associated trans clique also in vivo.

In all, we conclude that trans-C captures even weak and/or unstable yet biologically meaningful trans subnetworks associated with RNA biogenesis.

Loci strongly bound by DNA-binding proteins often exhibit significant interactions in trans

The identification of dense subnetworks of trans Hi-C contacts that are enriched for RBM20 targets supports our hypothesis that RNA biogenesis influences 3D chromatin organization by bringing into proximity coregulated nucleic acids, so as to maximize the efficacy and specificity of their processing (Fig. 5A; Bertero 2021). We specifically propose that, as in the case of RBM20, “RNA factories” arise from the clustered binding of trans-acting factors to one or more core coregulated genes and/or their encoded transcripts. These, in turn, recruit accessory targets of the same factors. This hypothesis predicts the existence of both transcription factories specialized for certain transcription factors (TFs) and/or chromatin regulators, and other RNA factories specialized for various RNA-binding and regulatory proteins. We set out to test this hypothesis systematically using trans-C, as an example of its potential applicability to address biological questions.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 5.

Trans-C identifies DNA-binding protein–associated trans cliques proximal to nuclear speckles. (A) Schematic of the mechanistic hypothesis for the formation of specialized RNA factories involving trans-interacting chromatin domains. Multiple copies of trans-acting regulatory factors (i.e., transcription or splicing factors) bind to core nucleic acids, aggregate to form new clusters and/or enrich pre-existing ones, and recruit accessory coregulated nucleic acids. RNA factories promote the efficacy and accuracy of RNA biogenesis processes (thicker black arrows). (B) Trans-C-identified subnetworks in lymphoblastoid cells built from loci characterized by strong binding of 110 DBPs have dense contacts. We plot the distribution of subnetwork weights for six types of sets of 40 loci: (1) loci with the highest number of ChIP-seq peaks for a given DBP (pink), (2) randomly drawn loci (gray), (3) top loci ranked by trans-C from a seed of five loci with the highest number of ChIP-seq peaks for a given DBP (orange), (4) top loci ranked by trans-C from a random seed of five loci whose starting subnetwork weight was matched to the seed of group 3 (yellow), (5) top loci ranked by trans-C from a seed of five randomly drawn loci (green), and (6) top loci ranked as for group 4 but starting from an interaction matrix that has been randomly shuffled (light blue). On average, sets seeded from loci most strongly bound by DBPs interact more strongly in trans than any of the other five types of sets of loci, including the stringent “matched random” control (P-values by Mann–Whitney U test). (C) For each DBP analyzed in B, we compare the weights of subnetworks obtained with trans-C from “top five DBP-bound” seeds (single data point) and “matched random” seeds (average of 1000 subnetworks ± SD). In red are comparisons with significantly different weights (P < 0.05 after FDR correction). Shaded areas highlight the top and bottom quartile of DBP-based subnetwork weights. (D) Visualization of selected significant DBP-associated trans cliques in lymphoblastoid cells, plotted as described for Figure 2D (PAX5 and MAX, strongest cliques; FOS and MLLT1, highest fold change of clique strength over average strength of cliques in the matched random null model). (E,F) Proximity to nuclear speckles of loci within trans-C-identified cliques, measured as the average SON and p-SC35 TSA-seq signal for the corresponding genomic regions. For each subnetwork, the signal is compared with that of an equal number of loci at the opposite end of the trans-C ranking. DBP-based cliques are overall significantly more proximal to both SON and p-SC35 than matched control sets (P-values by Mann–Whitney U test). Cliques that are significantly stronger compared with the null model in the analysis from panel C (Sign.) are in red, and nonsignificant ones (N.s.) are in blue. (G) The proximity to SON for significant versus nonsignificant cliques from panel C is significantly different by Mann–Whitney U test. (H) The strongest DBP-based subnetworks correspond to DBPs with a higher intrinsically disordered protein (IDP) score. We plot the IDP scores for DBPs resulting in the bottom and top quartiles of DBP-based subnetworks from panel C. The difference is statistically significant by Mann–Whitney U test.

First, we focused on DNA-binding proteins (DBPs), hypothesizing that the genes most strongly bound by a given DBP would be associated with strong TIDs. To test this notion, we used the most deeply sequenced Hi-C data set reported to date: an ultra-deep Hi-C map of human GM12878 lymphoblastoid cells (Harris et al. 2023). The ENCODE Project produced chromatin immunoprecipitation sequencing (ChIP-seq) data for 110 DBPs in this cell line (The ENCODE Project Consortium 2012), providing an ample resource to test our hypothesis in a systematic manner. For each ChIP-seq data set, we counted the number of peaks in each genomic bin.

First, we took the 40 bins with the most peaks for each DBP and calculated the weight of the subnetworks formed by these bins. The distribution of this subnetwork weight across all 110 DBPs is shown in pink in Figure 5B. For comparison, we randomly drew 1000 sets of 40 bins and plotted the distribution of their weight in gray. Clearly, the subnetworks of loci selected based on ChIP-seq peak density formed stronger interactions in trans than random sets of loci. This is a first important hint that many DBPs may be indeed involved in specific trans contact networks.

Second, for each DBP individually, we formed a seed by selecting the five bins with the most peaks from its corresponding ChIP-seq track, and we ran trans-C to identify a set of potential interactors in trans. We took the top 40 predicted bins for a given DBP and observed that these bins were enriched with ChIP-seq peaks not only for the DBP that spawned the seed but also for ChIP-seq peaks of other DBPs (Supplemental Fig. S4A). This is not surprising because many DBPs act in concert, and many loci contain proximal binding sites of several DBPs (Ibarra et al. 2020). When examining the weights of subnetworks formed by trans-C (“top five DBP-bound seeds,”) (Fig. 5B, orange), we noted that they were heavier on average than the subnetworks based on the ChIP-seq signal only (“top 40 DBP-bound loci”; pink). This observation validates that trans-C finds loci that interact even more strongly in trans with the seed bins than just the bins most bound by the respective DBP.

We also assessed how well trans-C can build dense subnetworks when it is seeded from biologically unrelated loci. To that end, we first drew 1000 times five random loci to use as seeds and ran trans-C (Fig. 5B, green). The subnetworks it built were significantly weaker than the DBP-based ones. This is likely because a randomly drawn seed set likely includes loci that are not interacting with one another, whereas the loci in the DBP-based seed tend to have strong interactions in trans. Thus, in order to establish a more stringent baseline, for each DBP we performed 1000 matched random controls with seeds of five random bins for each DBP (Fig. 5B, yellow). When using this matched random seed, the subnetworks that trans-C found were once again weaker on average than the ones it found using DBP-based seed (Mann–Whitney U test P = 0.006). At an individual level, the weight of subnetworks for 53% (58 out of 110) of the DBPs identified from the DBP-based seeds was significantly stronger than that from matched random seeds (P < 0.05) (Fig. 5C). Visualization of the two strongest subnetworks, associated to the TFs PAX5 and MAX, and the two most-significant subnetworks compared with their matched random controls, linked to the TF FOS and chromatin regulator MLLT1, demonstrated that these are interconnected with strong significance (Fig. 5D; Supplemental Fig. 4D).

Next, we performed an additional control in which matched random seeds were selected from a network that was previously randomly shuffled to remove all specific signals resulting from interchromosomal structure: As could be expected, the resulting cliques were on average the weakest recovered by trans-C (Fig. 5B, light blue), and individual comparisons for DBP-associated cliques were all significant compared with this type of control (Supplemental Fig. S4B). This confirms that interchromosomal genome architecture is far from random and that trans-C identifies signals much stronger than random noise. In all, these observations confirm the common sense conception that the interchromosomal interactions of biologically unrelated loci are mostly noise, while providing more rigorous support to the hypothesis that coregulated loci are often enriched for trans contacts (Bertero 2021).

To provide additional validation, we examined the strength of the DBP-associated cliques in an orthogonal data set based on split-pool recognition of interactions by tag extension (SPRITE) (Quinodoz et al. 2022), a proximity ligation-independent method to detect higher-order interactions within the nucleus. Specifically, we processed a SPRITE interaction matrix for GM12878 cells (Quinodoz et al. 2018) and evaluated whether the cliques’ trans-C identified using the Hi-C data exhibited strong interactions in this orthogonal SPRITE data set. All of the 110 DBP-based cliques showed much higher weight than a background of 1000 randomly drawn sets of loci (P = 6 × 10⁻⁴⁸, Mann–Whitney U test) (Supplemental Fig. S5A). To use a more stringent background model, we statistically assessed whether the trans-C-derived subnetworks for each of the analyzed DBPs were stronger than their “matched random” controls in the orthogonal SPRITE data (Supplemental Fig. S5B). We observed that the majority (30 of 58) of the DBP-based subnetworks that were significantly stronger than their matched random controls in the Hi-C data were also significantly stronger than their matched random controls in the SPRITE data (Supplemental Fig. S5C). Notably, our top five most-significant cliques were all significant in the SPRITE data. These findings suggest that trans-C reliably identifies cliques that exhibit strong interactions in trans also when these are measured by proximity ligation-independent sequencing approaches.

Lastly, we evaluated whether the loci that are part of a trans-C clique are physically closer to one another. To that end, we used orthogonal imaging measurements with multiplexed error-robust fluorescence in situ hybridization (MERFISH) (Su et al. 2020). We computed the trans-interaction proximity matrix for loci in the IMR-90 fibroblast cell line and downloaded the corresponding Hi-C matrix along with the ChIP-seq tracks for 16 DBPs available in the ENCODE portal. Because the MERFISH study involved only 1041 loci, we devised a twofold experiment. First, we subsetted the Hi-C matrix to the loci surveyed in the MERFISH study and ran trans-C using the five loci most bound by a given DBP as a seed. We calculated the average trans-proximity in the MERFISH data set for the loci in the subnetworks’ trans-C found using the Hi-C data and compared them to 1000 randomly selected sets of loci. We found that the trans-C cliques exhibited significantly higher proximity in the orthogonal imaging data set (P = 3 × 10⁻¹¹, Mann–Whitney U test) (Supplemental Fig. S6A). Second, we ran trans-C on the full Hi-C matrix with the five loci most bound by a given DBP as a seed to obtain a ranking of all loci, and then, we looked at the proximity of the 40 MERFISH loci that were ranked closest to the top of the trans-C ranking versus the proximity of the 40 MERFISH loci that were ranked closest to the bottom. We observed that for all DBPs the trans-C top-ranked MEFISH loci had significantly higher proximity than the bottom-ranked ones (P = 3 × 10⁻⁵, binomial test) (Supplemental Fig. S6B). These data extend the cross-validation of trans-C predictions of DBP-associated trans cliques from Hi-C data with orthogonal methods, including those based on imaging.

DNA-binding-protein-associated trans cliques are proximal to nuclear speckles

Our hypothesis is that DBP-associated cliques represent specialized RNA factories. To test this, we examined the nuclear localization of loci involved in trans-interacting subnetworks identified by trans-C. Specifically, we asked whether they are near nuclear speckles, membrane-less subnuclear organelles involved in various aspects of DNA and RNA metabolism (Galganski et al. 2017). To that end, we turned to data generated by tyramide signal amplification sequencing (TSA-seq) (Chen et al. 2018). TSA-seq is an experimental protocol that provides a “cytological ruler” for estimating mean chromosomal distances from nuclear landmarks genome-wide. We used the log₂ fold change of TSA-seq signal compared with an input control from the lymphoblastoid K562 cell line (Chen et al. 2018). This measurement captures the distance to a target protein from loci genome-wide, with higher values corresponding to shorter distances and lower values to longer distances. First, for each DBP subnetwork we calculated the average TSA-seq signal strength when probing the SON protein, which plays a crucial role in the formation of nuclear speckles. Indeed, the SON TSA-seq score is proportional to the cytological distance of genes from nuclear speckles, and it can be even calibrated to estimate mean distance in micrometers (Chen et al. 2018). Next, as a control, for each DBP we took the 40 bins ranked lowest by trans-C but containing at least one ChIP-seq peak, and we compared their average SON TSA-seq score to that of the trans-C-identified subnetwork (Fig. 5E). We observed a statistically significant shift to higher values for the trans-C subnetworks, supporting the conclusion that these loci are closer to nuclear speckles. An analogous analysis leveraging on TSA-seq data for phosphorylated SC35 (p-SC35), a splicing factor that also marks nuclear speckles, led to similar results (Fig. 5F). We noticed that cliques that were significantly stronger than their matched random controls appeared to be particularly close to nuclear speckle markers, particularly SON (Fig. 5E). Indeed, the SON TSA-seq score was significantly higher for these cliques compared with the other cliques (Fig. 5G). Collectively, these observations suggest that several DBP-associated cliques identified by trans-C represent RNA factories located at nuclear speckles.

Examining the distribution of the trans-C subnetwork weights identified for different DBPs (Fig. 5B, orange), we noticed a bimodal distribution, indicating that some DBPs are associated with stronger TIDs. This bimodality did not correlate with differences in the expression level of the two groups of DBPs or in their preference to bind to loci in the A or B compartments (Supplemental Fig. S4C). Intrinsically disordered regions (IDRs) within proteins, which lack a defined tertiary structure and are thus prone to self-aggregation, are emerging as an important mediator of subcellular condensates involved in multiple aspects of cell function, including nuclear regulations (Wright and Dyson 2015; Hirose et al. 2023). We thus investigated the correlation between the intrinsic disorder in DBP structure and the strength of the trans-C subnetworks they are associated with. For this analysis, we took the DBPs whose seeds gave rise to the strongest and weakest subnetworks (Fig. 5C, top and bottom 25% along the y-axis; Supplemental Table S4). We calculated the average intrinsically disordered protein (IDP) score for each DBP in the two groups (Mészáros et al. 2018), and plotted them in Figure 5H. The difference between the two groups was statistically significant (Mann–Whitney U test P = 1.4 × 10⁻⁵), suggesting that the DBPs with more IDRs form stronger interactions in trans. In all, trans-C allowed us to identify a large set of DBP enriched for IDR regions that are involved in strong TIDs proximal to nuclear speckles and that may thus be important in efficient transcriptional regulation of their target genes.

Selected RNA-binding proteins are associated to significant nuclear speckle–proximal trans cliques

Encouraged by the results on DBP subnetworks, we also examined whether RNA-binding proteins (RBPs) are generally associated with TIDs. Only a few RBP binding profiles are available for GM12878. Thus, for this analysis, we turned to K562 cells, another human lymphoblastic line, for which ENCODE reports 139 eCLIP data sets and a deep Hi-C matrix. Similar to our DBP analysis, we counted the number of peaks in each genomic bin for each RBP (mapping RNAs to the encoding DNA loci). Then, for each RBP individually, we formed a seed by selecting the five bins with the most peaks and ran trans-C to identify a set of potential interactors in trans. To form a null model per RBP, we repeatedly drew 1000 contact frequency–matched random seeds. We report the average total weight of the matched random seeds compared with the RBP seed in Figure 6A. Most RBP subnetworks built by trans-C were comparatively as dense as those from the corresponding matched random control, lying broadly along the y = x line. Nevertheless, several outliers were notably denser. To assess this observation quantitatively, we performed a signed ranked test per RBP and FDR controlled the corresponding P-values. Thirteen proteins had corrected P-values lower than 0.05, which we considered as a significance threshold (Fig. 6A, red; for list, see Supplemental Table S5). Distinctly from DBPs, RBPs associated with significantly stronger trans-C subnetworks were not characterized by higher IDP scores (Supplemental Fig. S7A), indicating that other characteristics may explain their specific behavior in trans genome organization.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 6.

Trans-C identifies RNA-binding protein (RBP)–associated trans cliques proximal to nuclear speckles. (A) A subset of trans-C-identified subnetworks in lymphoblastoid cells built from loci characterized by strong binding of RBPs have denser contacts than the corresponding matched random null model. We plot the weight of a RBP-based subnetwork and the average weight of 1000 matched random seed subnetworks (error bars correspond to the SD) for 139 RBPs. In red and listed by name are those with significant P-values after FDR correction. (B,C) RBP-associated cliques identified by trans-C are significantly closer to nuclear speckles (stronger SON TSA-seq signal) and significantly further away from the nuclear lamina (weaker Lamin A TSA-seq signal) than matched control sets at the opposite end of the trans-C rankings. (D) Visualization of selected significant RBP-associated trans cliques in lymphoblastoid cells, plotted as described for Figure 2D.

We repeated the analyses of TSA-seq data and observed a very strong and significant global correlation between the trans-C subnetworks and both SON and p-SC35 signal (Fig. 6B; Supplemental Fig. S7B). This indicates that the trans-C eCLIP subnetworks are close in cytological distance to nuclear speckles. This same trend manifested for POL1RE TSA-seq signal, a measurement of proximity to transcription factories (Supplemental Fig. S7C). We observed the opposite result when we examined the Lamin A TSA-seq signal, indicating that the trans-C subnetworks are significantly further away from the nuclear lamina (Fig. 6C). This finding is in line with the notion that sites of active RNA biogenesis are localized in the euchromatic nucleoplasm and away from heterochromatin regions associated with the nuclear lamina.

Visualization of significant RBP-associated subnetworks showed that these are noticeably interconnected (Fig. 6D), more so than DBP-associated subnetworks of similar strength (Fig. 5D); on the other hand, the individual pairwise interactions were less significant, possibly owing to a more transient nature. In all, trans-C allowed us to identify a subset of RBPs associated with nuclear speckle–proximal and transcription factory–proximal TIDs that may contribute to gene regulation.