2.1 Notations

In most methods, chromosomes are modeled as a serie of beads in 3D, each bead representing a loci or genomic window of given length. I denote by X∈R3×n the coordinate matrix of the 3D model, where n denotes the number of beads in the genome (for example, for P. falciparum at 20 kb resolution, n=1173). I denote by xi∈R3 the i-th bead’s coordinates, and by dij the euclidean distance between bead i and j

The contact map can be summarized by an n-by-n matrix C∈Rn×n, where each row and each column corresponds to a genomic locus, and each entry cij to the number of times loci i and j have been seen interacting. Note that the matrix is by construction square and symmetric (See Figure 3).

Figure 3:

Contact maps of P. falciparum’s chr7.

Both contact maps show enrichment of contact counts in VRSSM clusters (A. Ring stage [19] of IT strains, B. Trophozoite stage of 3D7 [20]) Dark blue corresponds to high interactions, light yellow to low interactions. The strong diagonal reflects the proximity of adjacent regions of the genome. The black line indicates the centromeres.

2.2 Consensus models

Consensus methods aim at inferring a model X∈R3×n that best represents the contact count matrix C∈Rn×n, usually through the minimization of an objective function O(X,C), sometimes under constraints.

The objective function O, sometimes called the “scoring function,” can be derived from known embedding algorithms (such as multidimensional-scaling methods), from statistical modeling of contact counts, or simply constructed in an ad hoc manner to fulfill a set of desirable properties.

2.2.3 Statistical models for contact counts

Another approach consists in modeling contact counts as random variables, with the 3D model being a latent variable. For example, Varoquaux et al. [23] propose an approach, called Pastis, that models contact counts as independent random Poisson variables where the intensity of the process between i and j is a function of the distance: cij∼Poisson(βdijα). The inference can thus be cast as maximizing the likelihood:

(2)maxα,β,Xℒ(X,α,β)=∑i <j⩽ncijαlogdij+cijlogβ−βdijα

Pastis can automatically adjust the parameters αβ of the counts-to-distance transfer function and infer a genome structure that best explains the observed data. The strength of this method comes from the robustness to low signal-to-noise ratio, thanks to a direct modeling of the noise through the statistical model (see Figure 4).

Figure 4:

3D models of P. falciparum’s genome architecture.

Using pastis-mds [23] (that implements a weighted MDS model) (A and B) and pastis-PO [23] (the Poisson modeling presented section 2.2.3) (C and D) on two Hi-C data sets of P. falciparum: a low signal-to-noise ratio data set at a Ring stage [19] (A and C) and a high signal-to-noise ratio data set at the Trophozoite stage [20] (B and D). pastis-MDS recovers a plausible structure for the high quality data set, but the non-uniform distribution of beads on low quality data set suggests that the method is not adequate on such a data set. In contrast, pastis-PO recovers plausible structures both on the high and low quality data sets. Each color corresponds to a chromosome. Large white beads mark centromeres, small blue beads telomeres and green beads VRSM clusters. Centromeric cluster is marked with a black line, telomeric cluster is marked with a dashed line.

2.3 Ensemble methods as a means to infer population of structures

Ensemble approaches aim at inferring a population of structures representative of the contact count map. The methods fall into two distinct categories: the first type casts the problem as a restraint-based optimization and samples local minima of the function [15, 22], whereas the second type proposes a statistical modeling of the problem and samples the posterior distribution [28, 31]. In short, the former is the ensemble version of MDS-based and ad hoc methods, while the latter is the ensemble version of statistical based methods.

2.4 Single-cell models

The last category of data-driven methods to infer the 3D architecture of the genome rely on a new protocol to probe single-cells for their 3D structures [37, 43]. Single-cell Hi-C is still in its early days, and, despite potential for assessing the variability of cell-to-cell genome architecture in a genome-wide fashion, only a handful of data sets are today publicly available. The contact maps originating from these data sets are very sparse, and specific methods to infer 3D structures need to be developed specifically for sc-HiC.

Akin to ensemble-local minima methods, the first approach is to consider each contact as a constraint and to formulate an under-constrained optimization problem. A population of structures satisfying the constraints can be found by sampling local minima [37]. Akin to consensus method, one can attempt to construct a distance matrix, either through manifold-based optimization [42] (by finding a low rank PSD approximation of the sparse contact map) or akin to Lesne et al. [29], by considering the weighted graph of interactions [38]. A classical MDS method applied to such a distance matrix then yields a consensus 3D model of the genome.

2.5 Model evaluation and comparison

A substantial difficulty in modeling the 3D structure of the genome is that model evaluation tends to be subjective. What is the relevant measure? “Truth” is generally not fully available, except for a few pairs of loci or in simulations. Is validating the colocalization of a pair or a few pairs of loci via FISH experiments enough? Are fit to the contact maps or agreements between modeling techniques relevant? Are the conclusions drawn from 3D models in agreement when these are inferred from different methods?

First, methods can be validated and compared against contact maps simulated from a known ground truth [23, 24, 29, 35]. Note that while this is a simple and natural first step for methods inferring a consensus structure, comparing and assessing robustness and accuracy of an ensemble of 3D models is much more challenging: in fact, it is still an untackled problem. Second, one can assess the stability and robustness of the inference with respect to (1) data bootstrapping; (2) contact map resolution (the models should not change as the resolution of the data varies) [23, 24]; (3) and in between biological replicates [23, 24]. As a cautionary remark, it is important to stress that it is not because a method is very stable that it is “good” with respect to any criterion: I can imagine a number of very stable methods, that would yet provide absolutely no insights on genome organization. Third, models can be compared to other sources of data, such as FISH [13, 20]. Fourth, biological plausibility of the resulting models can be considered: are the beads uniformly distributed in the cell [20]? Are known hallmarks of the genome architecture such as centromeres clustering preserved?

These are a handful of ways to assess plausibility and accuracy of 3D reconstruction, but many avenues in model evaluation and comparison are yet to be explored.

3.1 Building a yeast nucleus

The budding yeast S. cerivisiae’s 3D structure has been extensively studied, both through 3C-type studies [12, 13, 27] and through bio-imaging experiments [44]. The small size of its genome, the well-known hallmarks of its genome architecture and the availability of high resolution contact maps and FISH data set quickly led several teams to investigate the minimal set of constraints needed to reproduce the hallmarks of its genome architecture.

Tjong et al. [45] and Tokuda et al. [46] model S. cerevisiae’s chromosomes as flexible random fibers under a small set of constraints. While the exact modeling proposed by the three groups differ, the set of constraints can roughly be summarized as: (i) the chromosomes are constrained into a spherical ball representing the nucleus; (ii) centromeres are constrained into a spherical ball tethered to the nuclear membrane; (iii) telomeres are tethered to the nuclear membrane; (iv) rDNA is constrained into the nucleolus, represented as a spherical ball opposite to the centromeres, (v) a volume-exclusion constraint, preventing the fiber from occupying a space already occupied by the polymer. One can then simulate a large set of random structures fulfilling the constraints, and generate a “volume-exclusion contact map”, or “VE map”, from this population of structures, considering that two beads that are less then 45 nm apart in any of the structures form a contact. The Pearson correlation of the volume-exclusion contact map and the Hi-C one are highly correlated, demonstrating this small set of constraints fully explains the observed counts. In addition, the population of structures also explains FISH experiments previously published.

3.2 Building a P. falciparum nucleus?

For at least some stages, P. falciparum has a lot of genomic architectural features in common with the budding yeast: the centromeres are strongly co-localized at one end of the nucleus, telomeres are in physical contacts with one another, … Could these primary architectural features also arise from a population of constrained but otherwise random of structures? Adapting the set of constraints to match biological knowledge of P. falciparum, Ay et al. [20] showed that the resulting simulated VE map not only yielded lower correlations, but also failed to show the same features as the original contact count matrix. In particular, the VRSM genes which display domain-like enrichment in interactions do not appear in the simulated VE map (see Figure 5). Can we add constraints on VRSM gene clusters to improve correlations between true and generated contact maps? Adding a constraint on all beads considered as VRSM clusters is not enough: running 100 experiments yielded structures that did not fulfill all the constraints, demonstrating VRSM genes cannot all cluster together in a cell.

Figure 5:

Volume Exclusion Modeling.

Observed/expected matrices illustrate either depletion (blue) or enrichment (red) in contact counts for each pair of loci. A. Observed/expected map for volume exclusion modeling of P. falciparumB. Observed/expected map for Hi-C data.

4.1 Structure stability across time points, clustering and other variance analysis

The reader may well ask: “how sensitive are the resulting 3D models to initialization?” Taking the case of the P. falciparum, one may wonder whether structures from the same time points but from a different initialization are more alike than structures from different time points. A natural way to answer this question is to perform some dimensionality reduction technique, such as PCA and visualize whether structures from the same time points cluster with one another. Typically, features would then be the pairwise Euclidean distance of each structure, possibly subsampled to ease computation. Performing such an experiment on 1000 consensus structures inferred using the statistical model proposed by Varoquaux et al. [23], and available in the package pastis as pastis-PO, demonstrates that the results are more stable across initialization than time points (see Figure 6).

Figure 6:

Stability of structures across the life cycle.

PCA analysis of the population of structures obtained by running 1000 Pastis-PO on the three data sets of Ay et al. [20], corresponding to the Ring (Ay-Rin), Schizont (Ay-Sch), and Trophozoite (Ay-Trop) stages and the B15C2 data set of Lemieux et al. [19], at the ring stage (Lemieux-Rin) demonstrates that structures are more stable across initialization than across time points. Note that the Ring stages of Ay et al. [20] and Lemieux et al. [19] do not cluster, reflecting the centromeres strong colocalization in one of the data sets and not the other.

The second question is: “are models locally consistent?” A possible answer to this questions is to divide the structures into overlapping sub-structures ranging from 5 to 20 beads, and to compute the pairwise root mean squared deviation between segments across all structures [22]. Segments overlapping within a certain range for a large number of models can be assessed as locally consistent, while others should be labeled as highly variable.

The last question one can ask is: “how do the structures differ?” Tackling this question is very challenging, but can be reformulated to “are the hallmarks of interest conserved across structures of the same stage?” Ay et al. [20] and Lemieux et al. [19] both identified P. falciparum folded in very specific ways, with VRSM genes highly interacting. Ay et al. [20] also observed strong clustering of the centromeres, and enrichment in interaction at the telomeres. These observations can lead to a rigorous approach to identifying and quantifying whether families of loci cluster in the structures.

4.2 Chromatin compaction and chromosome entanglement

3D models can be used to estimate chromatin compaction and chromosome entanglement: distinguishing between open and close chromatin allows to relate the models to gene expression, open region being more accessible to the transcription machinery and thus genes more likely to be expressed. Chromatin compaction can be estimated by looking at the number of base pairs in a region defined either by volume if the scale of the structure is known, or by a percentage of the size of the structure if it is not [22]. Another idea is to sample random beads of a certain diameter, and assess how many loci are seen interacting. Applying this latter method on the three models of P. falciparum, Ay et al. [20] show that the trophozoite stage exhibits a more open chromatin than the ring and schizont stages. This finding is consistent with the transcriptionally active state of the parasite during this moment of the life cycle. A similar analysis, but counting the number of inter-chromosomal interactions, can help to assess the chromosome entanglement of the structure.

4.3 3D gene set enrichment

To assess whether groups of genes are colocalized in a 3D model, Ay et al. [20] leverage a statistical method developed by [47], which requires labeling each pair of loci in two groups: “close” or “far.” The authors used varying distance thresholds (10%, 20% and 40% of the nuclear diameter) to deem a locus pair “close” and labeled all remaining pairs in the set as “far.” The authors then compared the enrichment of loci pairs of a group being “close” and “far” by resampling loci among a same chromosome.

This approach dichotomizes loci pairs into two groups, and checks for the enrichment of a label in one of the two groups. Capurso and Segal [48] present an approach, called MPED, that avoids this step, and instead directly estimates the significance within the 3D model. Briefly, for a group G, MPED computes a test statistic:

where dij is the Euclidean distance between bead i and bead j. The null distribution is estimated empirically by resampling 105 times with preservation of the chromosome structure. If the M statistic is smaller than the mean of the null distribution, it is compared to the lower tail of the distribution and indicates co-localization. If the M statistic is larger than the mean of the null distribution, then it is compared to the upper tail of the distribution, and indicates dispersion.

Applying MEPD to the Trophozoite stage, Capurso and Segal [48] confirm that centromeres, telomeres, VRSM genes (both overall, subtelomeric and internal) colocalize.

4.4 Integrative analysis of gene expression and 3D structure using KernelCCA

Last but not least, an exciting contribution of Ay et al. [20] is the integrative analysis of gene expression and 3D structure using an unsupervised learning technique called “kernel Canonical Correlation Analysis” (kCCA) [49]. The goal of this analysis is to explore the relationship between gene expression and 3D structure, by extracting a set of gene expression components that exhibit coherence with respect to the 3D structure. While the components aren’t necessarily an actual gene expression profile, the genes of interest must somehow either be highly positively or negatively correlated with a component, and those genes should exhibit some form of coherence in terms of their 3D structure, like co-location. It can be helpful to think of this procedure as performing a principal component analysis (PCA) gene expression components extracted are correlated with the 3D structure.

Let us take a closer look at how kCCA is formally used in this context. Consider the set of n genes g∈G. Each gene g is represented on the one hand by its gene expression profile e(g)∈Rp and on the other hand by its 3D position x(g)∈R3. Assume the set of gene expression profiles is mean-centered and of unit variance.

The goal is to extract a gene expression component v∈Rp,∥v∥=1 such that it is both representative of set of gene expression profiles but also correlated with the 3D structure.

First, let’s tackle the question of representativeness of the set of gene expression profiles. We can identify a component v with such properties as to compute the percentage of variance explained by this component as follows:

Note that maximizing V(v) results in finding the first principal component v of the PCA. We can then define a score s∈Rn for each gene by computing the projection of each gene expression profile onto the component: s(g)=vTe(g). Any gene important to component v will by highly negatively or positively correlated with that component and thus either have a strongly negative or positive score s(g).

Now, let’s turn to the question of assessing coherence with respect to the 3D structure. Given a vector of scores f∈Rn, how can we assess the smoothness of these scores along the 3D structure? Ay et al. [20] leverage a standard approach in kernel methods, in which the smoothness of a score f is quantified by the function:

where K3D is the Gaussian kernel matrix of the genes’ 3D coordinates. The smaller S(f)is, and the more smoothly f is distributed in 3D.

So far, the two measures of representativeness and smoothness are independent from one another. However, if the scores s and f are required to be as correlated as possible, any genes highly correlated with each gene component v will also be co-localized in space.

To solve this problem, a common trick is to leverage reproducible kernel Hilbert space (RKHS) theory and cast the optimization in dual form. First, it can be shown that any candidate component v can be written as a linear combination of the gene expression profile: v=∑g∈Gα(g)e(g):α is called the dual coordinate of v. Let Kg∈Rn×n be the gram matrix of the gene expression profile, obtained by computing the inner product between all expression profiles: Kg(x,y)=e(x)Te(y). We can thus rewrite equation 3 as:

As K3D is invertible of dimension n×n, any score f can be written as f=∑g∈GK3Dβ, and the measure of smoothness as:

The optimization problem can now be cast in dual form as follows:

where λ is a penalization parameter. Solving this optimization problem identifies s and f such that; (1) the two scores are correlated with one another; (2) s maximizes V(s); and (3) f minimizes S(f). These optimization problems can then be solved efficiently using a generalized eigenvalue decomposition.

Ay et al. [20] apply this method using gene expression profiles and genes’ location extracted from the 3D models of the genome architecture and find gene expression profiles highly correlated with the 3D structure. Ranking the genes with their projections onto the gene components, Ay et al. [20] demonstrate that several gene families and Gene Ontology (GO) terms are enriched both close to the telomeres and at the opposite end of the nucleus. This method could easily be extended to other data set such as histone modifications, ATAC-seq, and beyond.