Global reorganisation of cis-regulatory units upon lineage commitment of human embryonic stem cells

Essential revisions:

Please address these main criticisms (more details are in the minor points section):

1) Many of the correlations between looping, gene expression and chromatin state are not novel. Focus on the new things, and present the data as a large dataset or resource for the community.

We thank the editors and reviewers for their constructive and thorough feedback on our study. We are happy for our work to be reconsidered as a Tools and Resources paper and we have taken on board the advice to focus more on the resource aspect of our study. We have achieved this, in part, by describing the large dataset more fully, adding further methodological details and support, and adopting a more cautious posture regarding the correlations between interactions and chromatin state.

2) The reviewers raised issues related to the statistical analysis and computational methods to determine PIRs (and their hierarchical clustering), TADs, and CRUs (below). Please address all of them.

3) Related to the previous point, the reviewers raised concerns about the claim that CRUs are distinct from INs, TADs. The main point to address is how the methods to determine these features compare and whether lack of overlap can be due to sensitivity of the computational methods.

We thank the editors and reviewers for their helpful feedback on these two points. We address these criticisms in our responses below.

Minor points:

1) One concern is that the ChICAGO method identifies significant interactions using only the general distance decay of interactions into account. It does not take into account that any pair of interactions between loci located within a TAD or IN, or between loci located in a similar compartment (A, B or subcompartment types) are generally higher than loci located in different TADs/INs/Compartments. As a result, a number of significant interactions are not necessarily "loops" (specific point-to-point interactions). This is not to say that these loci do not interact more frequently than expected given their genomic distance, but these interactions can reflect general higher order structures such as domains etc. This can explain why many interacting loci appear to not have chromatin marks: these can be "bystander" interactions that are the result of a nearby structural feature. I recommend that the authors repeat their analysis using a background model that takes domains into consideration.

The reviewers are correct to point out that CHiCAGO’s domain-unaware model somewhat overestimates the background read count (i.e. that expected at random) for TAD-crossing interactions (Author response image 1A). As a result of this, CHiCAGO is relatively more stringent (conservative) in calling TAD-crossing interactions than within-TAD interactions. In practice, the detected interactions are likely to be of comparable strength, as evidenced by the very similar read counts at the detected within-TAD and TAD-crossing interactions at a given interaction distance range (Author response image 1B). Therefore, the key observation of our study that TADs constrain, but do not fully insulate, promoter interactions is unaffected: the TAD-unaware CHiCAGO background only pushes us to err on the side of caution in detecting TAD-crossing interactions.

In the future, it would indeed be worthwhile to incorporate domain information into the CHiCAGO background model. However, this is a complex and challenging task. One of the challenges that need to be addressed, for example, is that the background depends not just on whether an interaction crosses a TAD boundary, but also on TAD boundary stringency, and does so in a way that scales non-linearly with distance (Author response image 1A). It is also likely that the background is affected by the number of crossed boundaries. Modelling all of these dependencies in the data must be performed with great care to avoid the danger of overfitting, which could lead to circular reasoning when attempting to make inferences about TADs.

Since the identified TAD-crossing interactions are at least no weaker than the within-TAD interactions, and since implementing TAD-aware interaction calling is likely a challenging project of its own, we feel that addressing this point in full goes outside the scope of this study.

Author response image 1

Download asset Open asset

Author response image 2

Download asset Open asset

2) The TAD calling procedure is quite simple. Is there any evidence that this procedure produces TAD calls of comparable quality to those produced by, for example, the HMM of Dixon et al. Nature 2012?

We agree that a comparison with Dixon et al., 2015 data is worthwhile and have now performed this analysis. This has revealed that the cutoff on HOMER boundary stringency scores that we initially used (∆Z>1) likely results in a more lenient TAD definition than the Dixon et al., 2015 HMM-based method, as evidenced by the higher numbers of detected TADs (4491 vs 3062) and smaller TAD sizes (Author response image 2A). Notably, the boundaries used in the initial submission showed a comparable enrichment for CTCF binding sites to Dixon et al., 2015 boundaries (Author response image 2B).

In the absence of a gold standard for TAD detection, it is difficult to judge which of the two methods is “more correct”. However, for consistency we now use a more stringent HOMER threshold (∆Z>2) that results in fewer TADs (2761 vs 4491 with ∆Z>1) and a comparable distribution of TAD lengths with the Dixon et al., 2015 TADs (Author response image 2A). More than 70% of our called TADs have at least one boundary that is consistently detected in the Dixon et al. study (Author response image 2C). To illustrate the boundaries called with both methods, we show their TAD calls for an example region in Author response image 2D.

The figures and the main text of the manuscript have been updated based on the new set of TAD calls. Using the more stringent set of TADs has resulted in fewer CRUs identified as crossing TAD boundaries (just over 20% vs ~40% in the initial set). However, this is primarily due to the larger sizes of the updated TADs rather than their higher ‘insulation power’. The ratio of observed-to-expected numbers of CRUs fully contained within TAD boundaries has in fact gone slightly down in the updated set (from 1.5-fold to 1.3-fold). We have also added the results for Dixon et al., 2015 TADs to Figure 3—figure supplement 1E, to enable a direct comparison. These data show that a similar proportion of CRUs cross Dixon et al., 2015 TAD boundaries (~27%) compared to those in our updated HOMER-defined set (~20%; ∆Z>2), demonstrating that the overall conclusion of this analysis is robust across at least two commonly used TAD-calling methods.

3) The TAD and IN comparison is unsatisfying. The main conclusion seems to be that CRUs are not particularly aligned with TADs and INs. This could be because the two phenomena are more or less independent of each other, or because the TADs, INs and CRUs are each defined in a fairly arbitrary fashion. If, for example, the simple TAD calling scheme employed here is not working well, then it would lead to the observed discordance.

We have now modified the TAD analysis as described above, which we believe addresses this point as far as TADs are concerned.

The IN calls have been taken verbatim from the study that introduced this concept (Ji et al., 2016). We noted, however, that the reported INs were often nested, so we have repeated the analysis using only the largest span of each overlapping set. As expected from the larger IN sizes, the numbers of IN-crossing CRUs went down somewhat, but still less than half of them (~45%) were fully contained within the extremities of IN boundaries (Figure 3—figure supplement 1I). Therefore, these results further strengthen a principal conclusion of our analysis that while INs constrain promoter interactions to some extent, they do not fully insulate them.

4) The finding that CRUs can extend beyond TAD or IN boundaries is interesting, but the data is not sufficient to claim a new "feature of genome architecture". Further analyses are needed such as: when a CRU extends beyond a TAD boundary, is this a weak TAD boundary? It is well known that TADs are nested structures. Do CRUs remain within the larger nested domains? This seems to be the case for POUF2. Quantification of boundary strength is required, as it seems not inconsistent with the TAD/IN models of enhancer action if the authors find interactions across weak boundaries, but not strong boundaries. Also, are these interactions beyond the boundaries equally strong, or are they much lower in contact frequency

To address this question, we plotted the distributions of HOMER strength scores for TAD boundaries that are crossed and not crossed by promoter interactions (Figure 3—figure supplement 1H). As can be seen, TAD boundaries crossed by interactions are indeed slightly weaker overall (Wilcoxon test p-value = 1.8e-14). However, the ranges fully overlap, and even the strongest TAD boundaries are not fully impenetrable to promoter interactions. We now also mention this analysis in the main text. Although these new findings have strengthened our analysis of CRUs, we take on board the reviewers’ concern about our interpretation of CRUs as new features of genome architecture, and have removed it from the manuscript.

As for the contact frequency of TAD-crossing and within-TAD interactions, it can be seen from Author response image 1B that it is very similar in both cases given the distance. In the short range, the identified TAD-crossing interactions even appear to have a slightly higher read count, most likely owing to the domain-unaware CHiCAGO background model, as discussed in our response to minor point 1.

Finally, we agree that empirically, TADs indeed represent nested structures, and tools for TAD calling in a hierarchical fashion are beginning to emerge (although neither the Dixon et al. method nor HOMER are among them). However, the conceptual disparity between the different levels of TAD “nesting” is still unclear. For example, from what level of the hierarchy should a sub-TAD be a considered a TAD, or a TAD be considered a “super-TAD”, etc.? It is possible to devise criteria that seem intuitively acceptable, but it will take time for the underlying assumptions of such criteria to acquire a solid structural and biophysical underpinning.

5) The main text claims that the promoter capture Hi-C data was processed to identify "significant cis-interactions." Looking at Methods, we learn that these contacts correspond to a threshold of 12, on a "log-transformed, weighted p-value" scale. This vague description makes it impossible to actually interpret the confidence associated with a threshold of 12, especially since we are told that the threshold corresponds to "a threshold of 5 in Chicago v1.0.0+ due to a soft-thresholding procedure introduced in this version." The choice of threshold must be fully justified and described. The text claims that Chicago uses an FDR control procedure, so the interactions that are presented here should be reported along with their FDR threshold.

We believe that, in part, this question has arisen because we made a confusing statement in the text. CHiCAGO’s weighted multiple testing treatment is not, strictly speaking, “a false discovery control procedure”, and we have now altered the main text to clarify this.

More generally, multiple testing correction when calling peaks in genomic data is a challenging problem, leading to the development of alternatives to standard approaches such as Bonferroni and FDR (such as, for example, IDR for ChIP-seq data). In PCHi-C data, an additional challenge is presented by undersampling, which leads to the violation of the underlying assumptions of both the IDR and conventional FDR approaches. In the case of FDR, this is because the majority of fragment pairs have only one or zero reads, leading to highly non-uniform p-value distributions. Furthermore, in PCHi-C data false discovery rates strongly vary with linear distance, partly due to the decreasing statistical power as distance increases. To account for this property, CHiCAGO adjusts the p-values using a weighting procedure, in which the weights trained on the reproducibility of interaction calls with distance (see Cairns et al., 2016 for details). Soft-thresholding was further introduced in CHiCAGO v1.0.0+ to shift the weighted -log-p-values (scores) such that the score of zero corresponds to the probability of a true interaction in the very short range given zero reads.

CHiCAGO scores are conceptually similar, albeit clearly not identical, to posterior probabilities resulting from an empirical Bayesian treatment. We view them primarily as a ranking measure, and use a heuristic approach based on integration with other types of genomic data to choose a threshold. Specifically, for a range of threshold settings, we compute the observed overlap of the identified PIRs with genomic regulatory features (such as specific histone marks) and the same overlap expected at random (estimated with random regions that are promoter distance-matched to the detected PIRs). In choosing a threshold, we balance the ‘unexpected’ overlap of PIRs with regulatory features with the total number of interactions detected at a given setting. As can be seen in Author response image 3A-B, at the selected threshold of 12 the relative difference in the overlap of PIRs with histone marks and CTCF versus random controls begins to saturate. Further increasing the threshold leads to a continued drop in the identified interactions (Author response image 3C-D) without a considerable further gain in the ‘unexpected’ overlap of PIRs with regulatory features.

We now explicitly mention the rationale for choosing the threshold in the Methods section. However, we agree that it is a crude and potentially subjective method. While varying the threshold within the reasonable range is unlikely to affect the principal conclusions of our study, we release the scores for each fragment pair as part of the Chicago objects on Open Science Framework, so the users can employ other methods of tuning the threshold, as well as apply non-threshold-based approaches to analysing our data.

Author response image 3

Download asset Open asset

6) The manuscript contains no details on precisely how ChromHMM was run. Details like software version numbers and any parameter settings should be given. Importantly, the text should specify whether the segmentation was carried out independently on each cell type, or jointly, and if the latter, whether the data tracks were "stacked" (i.e., 8 marks per position) or "concatenated" (4 marks per position, across two concatenated copies of the genome).

We apologise for this omission and have now added full details on the chromHMM procedure to the Methods section. Briefly, the segmentation was carried out jointly, and the data tracks were concatenated.

7) Details of the "curation" of the 16 inferred states should also be provided. When the 16 states were collapsed to four states, what was the basis for this collapsing? It seems strange that some states were eliminated at this stage because they contained multiple marks.

We have added full details on how the chromHMM states were curated to the Methods section. Briefly, the states were curated into “Active”, “Poised”, “Polycomb-associated” and “Background” based on the presence of active (H3K4me1/H3K27ac) and repressive histone (H3K27me3) marks. Overall, we were able to unambiguously assign one of these chromatin states to 81% of PIRs in both cell types. Two more curated states were defined, but not considered further:

1) “intermediate enhancers” characterised by H3K4me1 alone; the fraction of PIRs bearing this signature (~1%) was too small to analyse them as an individual category, and yet we felt that merging them with any other category was not biologically appropriate;

2) the ambiguous (“mixed”) state that showed enrichment for both H3K27ac and H3K27me3; we felt that this conflicting pattern likely arose from a technical issue such as heterogeneity within the samples.

8) How are PIRs defined: at the level of individual restriction fragments, or are adjacent fragments that both score significant merged into a single PIR?

PIRs are defined at the level of individual restriction fragments, with no signal smoothing or pooling across multiple sites. We have now clarified this in the Methods section. We have opted against further peak processing, because at present there is no reliable way to “fine-map” interaction signals without making a priori assumptions about the “causal” regions driving the signals, which we feel still require further mechanistic support at this stage.

9) The authors do not show many tracks of real data (only in one of the supplemental figures). Throughout they only show arcs to indicate significant contacts. This makes it difficult to assess the quality of the data.

We now provide two additional figures (Figure 1—figure supplement 4 and Figure 4—figure supplement 3), in which we plot the raw data for all examples shown in this study.

10) Promoters often interact with CTCF sites. Did the authors find this in their dataset as well?

We have now incorporated this analysis. Indeed, PIRs are highly enriched for CTCF binding sites, as we now show using ENCODE CTCF ChIP-seq data for hESCs in Figure 1—figure supplement 2B, and describe in subsection “Identification of putative regulatory elements and their associated gene promoters” of the main text.

11) The authors focus on cis-interactions. Can the authors describe trans interactions as well? Given that trans-regulation is considered rare, this analysis will provide context for interpretation of the cis data.

We now describe in the main text that we detect 338 trans-interactions in ESCs and 266 in NECs (representing <0.5% of total interactions). Consistent with our previous observations in other PCHi-C datasets (Cairns et al., 2016)(, a large proportion of trans-interactions (50.3% and 64.3%, respectively) reflect large-scale compartment signal that is detectable in conventional Hi-C data at 1-Mb resolution for the same cell types (data not shown). Given the small numbers of detected trans-interactions, we felt we should focus on cis-interactions in the manuscript. However, we agree that trans-interactions may be useful for future reanalyses of our data, and report them in the data files released on Open Science Framework.

12) Did the authors detect any known enhancer-promoter pairs? Also, the authors imply that the distal elements that touch promoters are regulating these promoters, but beside correlations between interactions, transcription and chromatin state of the distal elements, no evidence for direct regulation is provided. This would require deletion of the elements, e.g. by CRISPR. The Vista analysis of the POUF2 interactions suggests that many PIRs are enhancers. That is similar in information as the observation that PIRs overlap loci with histone modifications associated with enhancers. The key thing to validate is that these enhancers indeed regulate the interacting promoter. The fact that the pattern of activity of these enhancers resembles POU3F2 expression is encouraging but not definitive. The authors should either provide direct evidence for functional relationships or clearly discuss the extent to which the current data predicts such relationships even though the data currently is correlative.

This is an important point. We have added new text in the Discussion to highlight that our data set is predictive of regulatory relationships, but that extensive functional validation is required to infer direct regulation. In the absence of a comprehensive “VISTA-style” catalog of functionally validated enhancer-promoter relationships, we hope that large-scale resources such as ours will provide a strong impetus for future validation studies.

13) The authors examine the dynamic rewiring and recoloring of interactions as cells differentiate from ESCs to NECs. They report gain and loss of interactions and these are related to gain or loss of chromatin marks. This is interesting, but what needs to be tested is whether this is related compartment changes (see above). This is not to say these interactions are not relevant, but it would be important to determine how they relate to larger scale compartment changes. One interesting implication could be that in fact cell type-specific compartment changes are driven by altered interactions between functional elements that are active in a cell type-specific manner.

We have now done this analysis. While ~2% of 1-Mb regions transition between A and B compartments upon differentiation (in either direction), nearly all high-confidence rewired interactions detected in our data (>99%) map outside of these “dynamic” regions. We now mention this observation in the main text.

14) How do the authors interpret interactions that involve non-expressed genes?

This is an interesting point, which we now present more extensively in the Discussion. There are several possibilities to interpret interactions with non-expressed gene promoters. For example, it is known that at least in mouse ESCs, inactive genes and their regulatory regions form extensive topological networks (see, e.g., Schoenfelder et al., 2015b ). However, we do not intend to imply that all promoter interactions have a gene regulatory nature: some of them are potentially structural. Furthermore, additional and not yet fully understood mechanisms of topologically-mediated gene repression and priming may also be at play, particularly given the fact that a number of regions unmarked by conventional chromatin signatures currently emerge as having roles in gene regulation (Pradeepa et al., 2016; Rajagopal et al., 2016).

15) The hierarchical clustering of PIRs by the prevalence of different chromatin state labels is fairly uninformative. The fact that these different categories exist in the data is not surprising, and the claim that these eight labels somehow have "potential implications for understanding the logic of signal integration at promoters" is not well supported.

We agree that the majority of CRUs having consistent chromatin marks between PIRs and the promoter is relatively unsurprising. However, we feel that the classification we provide will be helpful for the users of our resource to make the most out of our data, and in particular to select targets for candidate-based studies. Additionally, our hierarchical clustering approach has been instrumental for detecting the relatively large numbers of ‘inconsistent’ CRUs, leading to the observation that the chromatin states at their promoters generally follow the ‘majority rule’ with respect to the chromatin states of the PIRs. We believe that these findings lead to specific and testable mechanistic hypotheses, and are therefore also interesting in this respect. Nevertheless, we concede that our claim about “potential implications[…]” was perhaps too far-fetched, so we have removed this sentence from the manuscript.

16) In the Abstract, the sentence "Here, we generate[…]" should mention what assay was employed.

This has now been done.

17) In the Introduction: "The extreme combinatorial complexity of Hi-C[…]" The complexity is really only quadratic, not combinatorial.

We have removed the words “extreme combinatorial” from this sentence to avoid confusion.

18) In the Results section: The URL for data availability should be included at the end of the sentence, "This data resource[…]"

We have now inserted a new sentence with details on the availability of raw and processed data.

19) Also in the Results section: "NEC PIRs were strongly enriched" This claim requires statistical support.

In the initial submission, the enrichment z-scores were given in Figure 2—figure supplement 1A, but not mentioned in the main text. We have now modified the text to rectify this.

20) New paragraph at "We next sought.[…]"

We have amended this.

21) Two of the observations offered in the Discussion do not seem to be supported by direct evidence in Results: "we find extensive promoter connectivity to regions associated with Polycomb-associated repression and poising, and "we.… detect large numbers of promoter interactions with regions devoid of chromatin features." These observations seem to be offered for the first time in Discussion.

Thank you for drawing our attention to this. We have now added details on these two points into the second section of Results (“Identification of putative regulatory elements and their associated gene promoters”).

22) The work does not show whether all the interactions are present simultaneously in a cell, or if the promoter state may vary between cells, depending on the probability with which the promoter is interacting with any of the potential alternative enhancers of a given chromatin state. Please discuss this issue in the main text.

This is a very good point. Indeed, Hi-C-based methods cannot distinguish between these two possibilities, and there is currently evidence in support of both the “hit-and-run” and the “chromatin hub” models. We discuss these topics in the section “Implications for developmental gene regulation by multiple enhancers”, which we have now expanded to explicitly outline these points.

23) Out of the ~21000 promoters, the authors focused their CRU analysis on only 16,000 protein-coding genes. Would the authors expect differences in the outcome if the other remaining promoters would be considered?

We now show the results for all captured promoters in Figures 1C and D, and have updated the main text accordingly. As expected, PIR enrichment for chromatin marks and the consistency between the promoter and PIR chromatin states are nearly identical when computed for protein-coding versus all captured promoters. Because of the level of similarity, we would also not expect any meaningful differences at the CRU level for all promoters compared to the currently analysed set of protein-coding genes. We would prefer to focus our analysis on protein-coding CRUs, however, as our mRNA-seq data interrogates only coding transcripts.

24) Out of the 16,000 protein-coding promoters about 9000 were defined as CRU. It may be good to know the selection criteria (i.e. is it based on a minimum number of PIRs per promoter?).

We apologise for not explaining this clearly. The ~9,000 figure includes all CRUs that have at least one PIR per promoter. We have updated the main text to explicitly state this fact.

25) Figure 3B defines a CRU as a set of PIRs that seem to be concurrently interacting with the promoter in the same cell. Maybe the authors could indicate in the caption that it cannot be ruled out that some PIRs may provide alternative rather then concurrent interactions.

We agree that this is a good idea and have amended the legend of Figure 3B accordingly.