The statistical-mechanics of chromosome conformation capture.

Since Jacob and Monod's characterization of the role of DNA elements in gene control, it has been recognized that the linear organization of genome structure is important for the regulation of gene transcription and hence the manifestation of phenotypes. Similarly, it has long been hypothesized that the spatial organization (in three dimensions evolving through time), as part of the epigenome, makes a significant contribution to the genotype-phenotype transition. Proximity ligation assays commonly known as chromosome conformation capture (3C) and 3C based methodologies (e.g., GCC, HiC and ChIA-Pet) are increasingly being incorporated into empirical studies to investigate the role that three-dimensional genome structure plays in the regulation of phenotype. The apparent simplicity of these methodologies-crosslink chromatin, digest, dilute, ligate, detect interactions-belies the complexity of the data and the considerations that should be taken into account to ensure the generation and accurate interpretation of reliable data. Here we discuss the probabilistic nature of these methodologies and how this contributes to their endogenous limitations.

Introduction

The fuzzy genome: From populations to single cells

The spatial organization of genomes has been investigated using different methodologies since studies into the cell’s components began. Techniques using low (e.g., refs. 1 and 2) and high through-put microscopy (e.g., ref. 3) and proximity-based ligation have enabled the generation of large data sets that facilitate the investigation of specific and global chromosome interaction patterns from an ever increasing number of organisms that includes bacteria,, baker’s yeast,, fission yeast, Drosophila, mouse, and human. Collectively these studies confirm that the positioning of loci varies within and between cells, even in identical conditions,, reviewed in refs. 13 and 14.

Proximity ligation assays are probabilistic in nature and as such rely upon population analyses. The probabilistic nature of these assays is due to both biological and technical aspects of the methods. Biological variation arises from the facts that: (1) DNA is highly flexible with a persistence length of 50 nm (150 base pairs). Therefore, a chromosome or DNA segment of tens or hundreds of megabases may assume many different folding conformations with equal probability and (2) despite the (unintentionally misleading) tendency to present proximity ligation data as unique conformations or summary heat maps, there is no a priori biological, evolutionary, or physical reason to assume a unique genome conformation in different individuals, or even within a population of cells in the same environmental condition(s). Technical variation also contributes to the probabilistic nature of the data, in this case variation arises from the cross-linking, digestion and ligation reactions which are inherently probabilistic. Thus, identifying an interaction does not mean the interaction is always present—rather that it is present within an undetermined proportion of the population. Conversely, the failure to detect an interaction does not mean that it never occurs, rather that the method does not detect it under those experimental conditions because the sites may not be cut or ligated.

Results from proximity ligation based experiments confirm contacts between genomic loci separated by large genomic distances (or between different chromosomes).-,-,- However, the existence of interactions does not automatically equate to biological significance, and interactions may simply originate from the very act of packing chromosomes into the nucleus or cell. Of course, this argument ignores the possibility that the restraints imposed by nuclear, or cell, shape have been integrated into chromosome organization over the course of evolution. Despite this, it is clear that some chromosomal interactions observed by modern techniques represent specific biologically relevant linkages,,,,.

Could there be a single spatial solution to a genome’s structure?

What is the average structure of a sportsperson playing a game? For the purpose of this analogy, any game will do; however we will specify that the game is rugby union. Over the course of a game, any one player’s structure changes as the game ebbs and flows toward its ultimate conclusion (Fig. 1A-E). The changes that occur to the player’s structure are not to the extent that arms or legs are removed—albeit sometimes they are broken—but rather the spatial arrangement differs from minute to minute according to his or her role within the current play. For example at the moment captured in Figure 1A the player’s (i.e., Ma’a Nonu) right hand is positioned next to his head as he participates in a ritualized challenge (e.g., haka; http://www.allblacks.com/index.cfm?layout=haka) performed at the beginning of the game. By contrast, at some point later in the game, the player’s right hand can be spatially adjacent to the left shoulder (Fig. 1B), in front of the right shoulder (Fig. 1C and D) or located at some distance (determined by the length of the arm) below the right shoulder (Fig. 1E). Could we determine the player’s structure by taking low resolution photographs of him or her at multiple time points across the game? Theoretically yes, the body parts would all be in similar relative positions (e.g., the head connected to the neck, arms to shoulders, etc.); however, even allowing for rotations that enable gross characteristics (e.g., the torso) to be aligned the average structure would be a dis-ordered amalgamation, the uniformity of which depends on the number of images which were taken and the body parts that were resolvable in the images.

Figure 1. Proximity ligation assays are probabilistic and determine the average structure within a population of cells. (A-E) Schematics of the individual conformations that Ma’a Nonu (All Black center) assumes during a game of rugby. Examples of the conformations Ma’a Nonu assumes are illustrated for: (A) the traditional All Black haka which occurs at the beginning of the game; (B) the act of passing the ball to his left.; (C) running; (D) a side-step and fend to avoid an opponent who is attempting to tackle him from the right; and (E) the act of scoring a try by diving to place the ball on the ground over the try line. (F) The average structural conformation obtained from (A-E).

This analogy represents the situation that would be predicted to occur for maps of genome organization that represent the different cell cycle stages (i.e., are unsynchronized) and even within cells that are synchronized to exactly the same moment of the cell cycle. Moreover, if we extrapolate the analogy to include the other players involved in the game, it demonstrates that averaging results across different cell types is also problematic.,

Variation is an inherent property of polymer structure

Precise positioning of entire chromosomes, which are in effect polymers consisting of several millions of monomers, would have too high an entropic cost to make it a feasible mechanism for controlling genomic processes in cells., In fact, it would imply holding the interphase chromatin together by some kind of rigid scaffold at support points not farther away from each other than the persistence length of chromatin—which has been estimated between 20 and 200 nm, a sizeable range but clearly much smaller than the extent of the typical eukaryotic nucleus. However, it is clear that the directed positioning of specific chromosomal regions does occur to different degrees during the cell cycle and in different cell types.-

Current experimental methods (reviewed in refs. 33-35) are limited and are unable to empirically define the precise positions of all loci within the genome of a single cell at high resolution. Therefore, to define features within the chromosome conformation ensemble, we are limited to modeling the ensembles of chromosome conformations using approaches developed for polymer physics coupled with data obtained from proximity ligation.,,,,, Despite the fact that there is inter-conformation variability between the structures in the ensembles, recent results have identified conserved patterns and clusters within the genomes that have been modeled to date.,,

Analyses of high resolution 3D structures of the S. cerevisiae genome that incorporated a proximity-ligation data set confirmed that yeast have preferred positions in the nucleus. Moreover, interaction-dependent clustering of tRNAs, early firing origins of replication and Gal4 upstream activating sequences were identified. Similarly, Ben-Elazar et al. reported the co-localization of co-regulated genes within space following analysis of statistically generated models that incorporated a global analysis of the S. cerevisise genome. These findings support the hypothesis that genome structure and function are interlinked.

Simulated three-dimensional (3D) structures of the Caulobacter crescentus genome, generated using 5C data and live-cell imaging, have illustrated that the C. crescentus genome is ellipsoidal with periodically arranged arms. Further clustering analysis of the individual structures within this ensemble identified four structurally similar configurations.

While the application of high-throughput sequencing to proximity-ligation data has led to significant increases in data generation and a leap in our understanding of genome organization, the reconstruction of 3D genome conformations based on such data is in its infancy. Appropriate methods for analysis are still being developed. The main problems in the linking of data to models of genome organization are that: (1) the ligation probabilities are not exclusively a function of the proximity, but influenced by many other factors such as intrinsic reactivity, local coverage by proteins, crowding, etc.; (2) these methodologies have the intrinsic property that they only provide information on pairwise interactions–simultaneous interactions of several loci in a cluster can only be inferred, not be directly identified; (3) restriction enzyme choice affects the frequency of detection of restriction fragments and ligation products; (4) PCR biases affect the chances of detecting and identifying interactions; and (5) the copy number of genomic regions affects the chances of detecting interactions with and between these loci. Despite these methodological limitations and the structural variability that exists within conformation ensembles, the results of 3D modeling have identified conserved patterns and clusters within those genomes that have been analyzed thus far. These models have the potential to generate new hypotheses about the relationship between nuclear structure and function.

Proximity ligation assays,,, have demonstrated that, in addition to the directly adjacent neighbors, loci interact with up to 38 partners depending on the organism or data set under study (Table 1). The variation in observed numbers of interacting partners for loci within an organism reflect the: (1) physical characteristics of the chromosomes; (2) depth of interaction coverage: e.g., the depth of sequencing, the choice of restriction enzyme, and resulting fragment lengths affect if an interaction is detected and hence the coverage; (3) relative complexity of the functional restraints that result in the formation of an interaction; and (4) the environment within which the genome was analyzed. From a functional perspective, it is tempting to speculate that an association with large numbers of partners indicates that loci are clustered at hubs, for repair, replication or transcription.- However, the validity of this interpretation is debatable as it ignores biological issues including the fact that cell volumes, and hence global chromatin compaction levels, vary and this variation can alter our ability to detect interactions by unblocking restriction sites. Moreover, the interpretation also ignores the technical issue that proximity ligation assays determine the population average and do not identify the combination of interactions that occur within a single cell.

Table 1. Loci that interact within a genome tend to do so with more than one locus. No. of interactions per fragment, the number of interactions that were empirically determined for fragments that interact with at least one other fragment within the data sets

			No. of interactions per fragment
Organism	Condition/cell line	Method	Mean	Max	Reference
Saccharomyces cerevisiae	Glucose	GCC	11a	30	17 , 21
	Glycerol lactate	GCC	6a	20	17
	Galactose	GCC	3a	21	17
Escherichia coli	exponential	GCC	7.83	38	5
	serine hydroxamate	GCC	6.67	25	5
Homo sapiens	GM12878	5C	3.88b3.25c	20	15
	K562	5C	4.07b3.84c	> 9	15
	Hela-S3	5C	5.35b4.62c	> 9	15

aMode number of interacting fragments for restriction fragments that interacted with more than just adjacent fragments within the S. cerevisiae genome. bMean number of interacting fragments for transcription start sites (TSS), throughout the ENCODE pilot regions representing 1% of the human genome, with at least one non-adjacent interaction for expressed genes in GM12878. cMean number of interacting fragments for TSS sites for non-expressed genes in GM12878. GCC, genome conformation capture; 5C, 3C carbon-copy.

At present, most of the existing proximity-ligation data are a population average of both large numbers of cells and and cell cycle stages. It could be argued that proximity ligation assays performed on single cells would enable the direct determination and identification of clusters within each cell that was being tested. However, this approach would still fail to clarify the problem completely for the following reasons: (1) restriction digestion is incomplete; (2) ligation is a probabilistic event that covalently joins only two free DNA ends; and (3) there is PCR and copy number bias.

Restriction digestion is a source of variation

The choice of enzyme affects the patterns that are observed and the interpretations that are made (e.g., ref. 5 vs. ref. 44). This is particularly important given the recent trends to map transcriptional activity back onto interaction networks, which requires small fragment sizes to allow the accurate overlapping of interactions with transcriptional or mediating factors (e.g., ref. 5 vs. ref. 44).

Recent experiments have shown that the choice of the restriction enzyme alters the amount of chromatin that is solubilized, thus affecting the frequency with which interactions are detected and refining the interaction pattern. For example, the ability of a fragment to circularize and be solubilized, is dependent on the length of the DNA fragment and hence the choice of the restriction enzyme. Moreover, restriction enzyme choice affects ligation rates with blunt end ligation is less efficient than sticky end ligation.-

The specificity of the enzyme and any star-activity needs to be taken into account, particularly given the large quantities of enzyme used in the standard protocols. Non-specific cleavage and subsequent ligation can complicate analyses. Finally, the efficiency of chromatin cleavage is important, as too low an efficiency will result in the detection of only adjacent interactions due to failure to separate the restriction fragments. Significant locus specific differences do exist, for example cleavage efficiencies of 85% and >70% can be achieved with HindIII and MboI within the β-globin locus. Similarly, cleavage efficiencies of between 65–89% are achieved using AseI at the mouse immunoglobulin kappa locus.,

Variable ligation is a significant cause of methodological variation in proximity ligation assays

Intra-molecular ligation is the critical step in all proximity ligation assays. As such, it is important to incorporate the physico-chemical consequences of this step into the interpretation of results from such assays. The process of ligation is an inherently probabilistic event. The probabilistic nature of this process results from the fact that ligation events mediated by enzymes, for example the T4 ligase, depend upon two compatible free ends simultaneously associating with the active site of the ligase. The association of all three components (i.e., two free ends and a ligase molecule) is the limiting event, hence the improvement in ligation frequencies achieved in reactions containing a 3:1 molar excess of insert to vector DNA during cloning, and enzymes that can bind and not release the free ends.

In the case of the intra-molecular ligations detected by proximity ligation, the local concentration of ligation competent free ends is artificially increased through the formation of intermolecular cross-links that hold the interacting restriction fragments together. However, the efficiency of the resulting ligation reactions is also affected by the relative flexibility of the DNA fragments. The flexibility of DNA for fragment lengths above 200 bp, is well known and given by the persistence length of 50 nm. However, the fact that the restriction fragments are (1) different lengths (Fig. 2), (2) bound by different proteins or protein complexes (e.g., nucleosomes or nucleoid associated proteins) that constrain their flexibility and ability to interact and (3) that the positions of the bound complexes or inter-linking are unlikely to be evenly distributed along or on the restriction fragments meaning that the length of DNA which is free to bend at the end can be limiting (Fig. 2B-D). There is ongoing debate as to the minimum length that can circularize.,,

Figure 2. Restriction fragments lengths and protein binding can limit the ability of free ends present on cross-linked fragments to ligate. (A) The circularization of linear DNA fragments is length dependent with fragments below ~70 bp classically regarded as not being able to circularize. (B-D) Cartoons illustrate how nucleosome (gray disc) binding at or near the end of restriction fragments can prevent the ligation of free ends and limit the potential ligation products from two restriction fragments (red and black lines).

Arguably the single most important limitation of the proximity ligation step is that any one restriction fragment can only ligate to two other restriction fragments—one at each free end. The chances of this occurring are very low, for example Sexton et al. only identified 2 products in 167 as consisting of three ligated restriction fragments that definitely originated from a tripartite complex. Moreover, even if only two ligation-competent fragments are cross-linked, ten different products can form during the ligation reaction (Fig. 3A). In fact the number of individual products that can form through ligation events within clusters containing incrementing numbers of restriction fragments (n) increases super exponentially such that complexes containing 10 ligation competent partners can produce 23 758 664 096 products (Table 2; Supplemental Material). Critically, not all of the products that form are able to be distinguished from each other, particularly when they are fragmented for sequencing (Fig. 3B).

Figure 3. Ligation of restriction fragments within a cross-linked complex results in super exponentially incrementing numbers of ligation products. (A) Ligation events between two restriction fragments, held together within a cross-linked complex (gray oval), results in 10 possible combinations of products. (B) Ligation events between three or more restriction fragments result in overlapping sequences following fragmentation for sequencing.

Table 2. The number of ligation configurations that are possible when 2–10 ligation competent fragments are held within one cross-linked complex

Number of interacting fragments	Number of different ligation configurations
2	10
3	76
4	764
5	9496
6	140 152
7	2 390 480
8	46 206 736
9	997 313 824
10	23 758 664 096

Recent experiments have shown that those fragments which do interact typically do so with more than just their adjacent fragments and up to ~30 partners have been identified within studies in populations of yeast cells (Table 1). As stated earlier, proximity ligation assays do not allow the identification of clusters within a single cell. However, measurements based on transcriptional activity and numbers of active polymerases have led to estimates that transcription factories contain eight polymerases engaged on different loci in Hela nuclei. Moreover, the formation of centromeric and telomere clusters indicates that biologically relevant clusters do exist. This complicates attempts to calculate actual ligation frequencies, and hence inferences about the percentage of the population in which an interaction occurs.

Repetitive elements also pose a potential bias in analyses of proximity ligation data as they are difficult to position on the genome with any certainty. Therefore, it is difficult to determine which particular repeat(s) are involved in a specific interaction. For example, telomeres and sub-telomeric regions are repetitive, have strand specific base compositions and are known to cluster in a wide a variety of organisms and at specific cell cycle stages (reviewed in ref. 60). This raises issues around the modeling of regions that are effectively invisible to the proximity ligation methodology. However, the junctions of specific repetitive regions can be mapped and, assuming they interact with other regions, can be used to localize the beginnings and ends of the repetitive domain. The structure of the remainder of the repetitive domain can be approximated using the biophysical characteristics of the DNA. However, the failure to be able to uniquely position interactions with repetitive regions is a limitation that needs to be borne in mind when analyzing and interpreting proximity ligation data.

Interactions with repetitive elements that are arranged in tandem (e.g., the Saccharomyces cerevisiae rDNA repeats) or episomes (e.g., mitochondrial genomes or high copy number plasmids) can be collapsed to one particular region or DNA element. This is a common approach, as current analyses of diploid and replicated genomes are unable to discern which copy of a particular chromosome is involved in an interaction—yet there is no a priori reason to assume that both copies are involved in the same interactions.

Current proximity ligation methods are practically and technically incapable of identifying all of the products that form from the complexes that are present within the individual cells. Therefore, the full complement of interactions that are occurring within a single cell cannot be identified even if all of the cells within a set population are individually assayed because the ligation products that represent some of these interactions cannot form—physically or statistically.

Bias due to PCR amplification

The amplification of sequences prior to sequencing, either as part of the proximity ligation methodology itself or as part of the library preparation for sequencing, can introduce bias into the results. This is particularly true for samples that have high nucleotide base biases (e.g., skewed AT or GC compositions,) and is dependent upon the enzymes, processes and equipment that are used to prepare the samples for sequencing. For example, bias can be introduced at any step during Illumina sequencing, including: (1) library preparation which can introduce intermolecular ligation events, (2) cluster amplification which can affect GC bias, (3) synthesis, and (4) post-sequencing data processing.

PCR bias affects both the chances of seeing a sequence, the sequencing error probability, and thus the chances of mapping it back to the reference genome., Protocols are constantly being modified to reduce these biases (e.g., ref. 61 and “PCR-Free” sample preparation kits [TrueSeq DNA PCR-Free Sample preparation kit; FC-121–3001; Illumina]). These alterations to sequencing protocols introduce differences that affect cross-study comparisons between existing and new data sets, or for data that is generated on different platforms at different centers. Such biases are obvious in run-specific frequencies for the identification of intermolecular ligation rates such that different numbers of intermolecular ligation events are observed in different preparations of the same sample.

It could be argued that the AT or GC biases we are discussing are not important when looking at genome spatial organization. However, not only do different genomes have different GC biases, but repetitive elements and regions with skewed AT or GC compositions exist within prokaryotic, and eukaryotic genomes, are biologically important and may play significant roles in genome organization. Centromeres are a classic example, however long nucleosome-free regions (LNFRs) in resting human T cells have skewed AT (average 73% AT) or GC (average 76% GC) contents, and these regions have been suggested to be involved in global remodeling.

With the exception of the terminal fragments on linear chromosomes, all loci interact with two linearly adjacent restriction fragments when assayed by proximity ligation. These linearly adjacent interactions complicate interpretations as they may represent: (1) real interactions mediated by some protein complex; (2) protection of the restriction site; or (3) incomplete digestion. As such, they are typically removed from subsequent analyses. While effective, the question remains whether the bioinformatic removal of these interactions from data sets overcomes all of the complications associated with their presence. Empirical methods, in particular the use of modified blocker primers, can be implemented to prevent or reduce the amplification of the adjacent restriction fragments in 4C assays and other few against many targeted proximity ligation assays. These methods result in the relative amplification of the signal from the less common long-range or non-adjacent interactions.

Bias due to copy number variation

Copy number of fragments or ploidy within the genome is not just a complicating factor for repetitive regions, rDNA repeats, plasmids or organelle genomes. Rather it is an issue for interactions involving any loci or chromosomes that exhibit copy number variation: (1) as result of genome replication; and (2) due to specific variation in response to environmental cues (reviewed in ref. 68). Variation due to replication can complicate analyses of populations with altered proportions in the different cell cycle stages. Alternatively, replication complicates bacterial analyses due to the fact that these cells show varying states of aneuploidy due to the existence of multiple instances of initiation from their origin of replication (e.g., Escherichia coli). As a result of these issues, copy number variation should not simply be considered as an artifact to be removed because it complicates the mapping of the connected loci. Rather, it remains possible that connections between loci that show copy number variation are in fact biologically significant and informative.

The ability to distinguish or not effects that are correlated with alterations to copy number is an acute issue for proximity ligation studies which do not incorporate some aspect of copy number variation analysis. The GCC technique enables copy number across the entire genome to be accurately determined and accounted for in studies of genome organization. While this is useful, the corollary is that the amount of sequencing required is currently prohibitive for studies of large genomes (>300 Mb). It may be possible to account for copy number in other global proximity ligation techniques using external standards that are applied at known concentrations prior to the steps that are used to enrich for ligated fragments prior to sequencing.

It is possible to determine the number of different configurations of interactions using an advanced inductive combinatorial argument (Supplemental Material). The outcome is an algorithm which gives the number of configurations a(n) as a function of the number of candidate fragments n, which can join from either end.

Conclusion

Given the intrinsic weaknesses within proximity-ligation assays, it is clear that the amount of information that can be obtained by direct analyses of identified interactions is limited. Moreover, while some of the biases can be compensated or corrected the direct interpretations of the crosslinking data directly will only inform on those elements which are physically connected and thus cross-linkable. It remains debatable how to proceed to overcome this limitation. Is it simply a question of depth of sequencing of the ligation products, or does the answer lie in the incorporation and analysis of in silico reconstructions or forward predictive models. Reconstructions of varying complexity have already been generated,,,,, and the analyses of these models is complex. Predictive models can be used to ask different questions and could be generated beginning with ensembles of polymer models that recapitulate the random unpacking of DNA within the nucleus., The integration of predicted chromosomal contacts generated using transcriptional, replication and repair networks would morph these random models into “accurate” testable representations of in vivo genome organization. Early attempts at this approach have been made., However, the combinations of paramaters that need to be used to generate “accurate” predictive models remain unknown.

Ambiguity comes not just from the visual representations of genomic organization; rather, it is a fundamental aspect of the proximity based ligation methodologies that are being used to study structure and the spatial organization of genomes themselves. Without methodological improvements that include the continued integration of alternative data, we will remain unable to accurately identify the pattern of structural spatial associations within single cells. Statistical limitations mean that the required improvements are unlikely to be based on the proximity-based ligation methodologies as they are currently employed. However, it remains undeniable that the potential advances in our understanding of the structure: function relationships within the nucleus are key to understanding not only nuclear processes but also the processes of development and how cells respond to their environment.