Simulation-based inference of yeast centromeres

TL;DR

This paper introduces a novel stochastic method for inferring yeast centromere locations using Hi-C data and simulated contact maps, addressing the challenge of unknown centromere positions in many species.

Contribution

It presents a new approach combining experimental Hi-C data with simulations to accurately locate centromeres in yeast genomes.

Findings

Successfully infers centromere locations in budding yeast

Combines experimental and simulated data for improved accuracy

Addresses limitations of previous point-estimate methods

Abstract

The chromatin folding and the spatial arrangement of chromosomes in the cell play a crucial role in DNA replication and genes expression. An improper chromatin folding could lead to malfunctions and, over time, diseases. For eukaryotes, centromeres are essential for proper chromosome segregation and folding. Despite extensive research using de novo sequencing of genomes and annotation analysis, centromere locations in yeasts remain difficult to infer and are still unknown in most species. Recently, genome-wide chromosome conformation capture coupled with next-generation sequencing (Hi-C) has become one of the leading methods to investigate chromosome structures. Some recent studies have used Hi-C data to give a point estimate of each centromere, but those approaches highly rely on a good pre-localization. Here, we present a novel approach that infers in a stochastic manner the locations of all centromeres in budding yeast based on both the experimental Hi-C map and simulated contact maps.

Figures

Equations

$$\hat{\theta }=\underset{\theta \in \Omega }{\mathrm{argmax}}\log p(C_{\text{ref}}|\theta ).$$

$$\mathbb{E}\left[\theta \mid \Vert S(C)-S(C_{\text{ref}})\Vert \le ϵ\right]\underset{ϵ\to 0}{\to }\mathbb{E}\left[\theta |C_{\text{ref}}\right],$$

$$\mathbb{E}\left[\theta |C\right]=\underset{S\in \mathcal{F}}{\mathrm{argmin}}\mathbb{E}\left[\Vert S(C)-\theta {\Vert }_2^2\right],$$

$${\hat{\mathcal{L}}}_{\text{DNN}}(\varphi )=\frac{1}{N}\sum _{1\le n\le N}\Vert S_{\varphi }(C^n)-{\theta }^n{\Vert }_2^2.$$

$${\mathbb{E}}_C\left[D_{\mathrm{KL}}\right(p(\cdot |C)\Vert p_{\psi }(\cdot |C)\left)\right].$$

$${\hat{\mathcal{L}}}_{\text{NPE}}(\psi )=-\frac{1}{N}\sum _n\log (p_{\psi }({\theta }^n|C^n))\text{ , }({\theta }^n,C^n)\sim p(\theta )p(C|\theta ).$$

$${\theta }^n\leftarrow {\bar{\theta }}^k+ϵ\text{ with }ϵ\sim \mathcal{N}(0,{\sigma }^2\text{Id})\text{ for }k=n\text{ mod }M\text{ and }n=1,...,N$$

$$w^{m,t}=\frac{\pi ({\theta }^{m,t})}{{\sum }_{k=1}^M{w}^{k,t-1}K({\theta }^{m,t};{\theta }^{k,t-1})}$$

$${\hat{\mathcal{L}}}_{\text{DNN}}(\varphi )=\frac{1}{N}\sum _{1\le n\le N}\Vert S_{\varphi }(C^n)-{\theta }^n{\Vert }_2^2$$

$$p(\theta \mid \Vert S_{\varphi }(C)-S_{\varphi }(C_{\text{ref}})\Vert \le ϵ)\to p(\theta |S_{\varphi }(C_{\text{ref}}))\approx p(\theta |C_{\text{ref}}).$$

$${\hat{\mathcal{L}}}_{\text{NPE}}(\psi )=-\frac{1}{N}\sum _{1\le n\le N}\log p_{\psi }({\theta }^n|S_{\varphi }(C^n))$$

Full text

1 Introduction

Hi-C maps have become one of the main assets in understanding DNA folding, notably through the study of chromatin loops or topologically associated domains (TADs) in mammalian cells tads ; loops . These maps capture the contact counts between fragments of chromosomes among a population of DNA into a 2D squared and symmetric matrix made of cis- and trans- blocks of interactions.

Beside chromatin loops or TADs, centromeres are also of great interest to genome structure investigation due to their essential role in many biological processes: they facilitate chromosome segregation through the formation of the kinetochore kinetochore during mitosis and meiosis, and act as key regulators of genome stability via the prevention of chromosome breakage. In yeasts, centromeres are highly compact regions spanning about 125 base pairs (bp) centrosize and tend to cluster near the spindle pole body within the nucleus. This clustering results in a distinct peak in the trans-contact counts Hi-C matrices, centered at the position of each centromere pair. Many studies have attempted to annotate yeast centromeres, usually through Fluorescent In Situ Hybridization (FISH) FISH or chromatin immunoprecipitation (ChIP) CHIP . However, these approaches often remain imprecise and may even fail to infer centromeres for some species inferfail . To bypass these limitations, newer methods have been proposed, which use Hi-C contact maps nelly ; nelle . Working directly with a Hi-C contact map or with the corresponding Pearson correlation matrix, they fit a Gaussian to each interaction peak to precisely infer centromere locations. Beyond the optimization of a non-convex function, these methods highly rely on good pre-localization of the centromeres to be precise and output only a point estimate of each centromere whereas it is actually a whole segment of chromosome.

We propose a novel approach to infer centromere positions that differs from existing methods in two key aspects. Firstly, we adopt a stochastic approach by quantifying the uncertainties about the centromere candidates we infer. Secondly, the inference processes are not only based on a reference Hi-C matrix (denoted $C_{\text{ref}}$) but also on simulated contacts maps (denoted $C$). We adopt a Bayesian approach where centromere positions (denoted $\theta$) are sampled from a prior distribution and the contact maps $C$ are generated from a custom-designed simplified simulator. We estimate the posterior distribution $p(\theta|C_{\text{ref}})$, which amounts to solving the following inverse problem:

Given contact map $C_{\mathrm{ref}}$, what are the most probable centromere positions $\theta$ to have generated it?

2 Methods

2.1 Framework

Setting.

We work with the budding yeast Saccharomyces cerevisiae for which the positions of its 16 centromeres and length of each chromosome (in base pairs) are known. We used data from Duan et al. duan to construct the reference Hi-C contact map. Centromere candidates $\theta=(\theta_{1},\ldots,\theta_{16})$ are sampled from a prior distribution $p(\theta)$ that is poorly informative, e.g. a multivariate uniform distribution where the interval is the length of each chromosome in each dimension. To simulate contact maps, we designed a simplified simulator (described below) that takes $\theta$ as input and directly outputs realistic $C$ without simulating any DNA folding. Other studies have used simulators to do Bayesian inference of biological elements arbona , but those biological simulators try to simulate the 3D folding of the chromatin in the nucleus before computing the corresponding contact map. This renders them too slow and unnecessarily complex if we want to have many contact maps in a reasonable time.

Contact maps and data normalization.

The contact map $C$ projects the information contained in a population of 3D chromatin foldings into a 2D squared and symmetric matrix made of cis- (or intra-chromosomal) and trans- (or inter-chromosomal) blocks of interactions between pairs of chromosomes. To construct it, we cut each chromosome into genomic windows of a given length (called resolution), e.g. 32 kilobases (kb). Each matrix entry is then a non-negative number, called the contact count, representing the number of times a given window was in contact with another one over the population (see Appendix A and Figure 3 for more details).

During inference, we use a reference Hi-C map $C_{\text{ref}}$ and simulate synthetic contact maps $C$. Hi-C contact maps have many biases due to sequencing and mapping errors or to the inherent structure of the chromatin hicnorm . Therefore, $C_{\text{ref}}$ is actually a normalized Hi-C map, where the normalization corrects those biases, iteratively forcing all rows and columns to sum up to one hicnorm . The quality of the contact map depends on the chosen resolution and the signal-to-noise ratio gets smaller if we work at higher resolution. We thus choose to set the contact maps at resolution 32 kb. In yeasts, the main informative part about centromeres rely on the upper trans-contact blocks: the matrix of contacts between chromosomes $i$ and $j$ contains an enrichment of interactions at the location of both centromeres $(\theta_{i},\theta_{j})$.

Simulator.

We exploit the structure of yeast contact maps to design a very efficient simulator that directly creates the upper trans-contact blocks given its centromeres positions $\theta$. At the centromere positions, the chromatin has a brush-like organization: chromosomal regions near the centromeres often enter in contact over the population whereas the further we move away from the centromeres, the rarer the contacts become. To mimic this effect, we simulate a Gaussian spot at the position ($\theta_{i},\theta_{j}$) for each trans-contact block. Between chromosomes, we also observe rare interactions over the population that we reproduce by adding Gaussian noise to all the trans-contacts blocks up to $10\%$ of the maximal contact count (see Appendix B with Algorithm 1 and Figure 4 for more details).

2.2 Simulation-based inference

Our goal is to infer $\theta$ from $C_{\text{ref}}$ using a probabilistic framework based on simulations. The usual way for doing so would be to search for the most appropriate $\theta$ for a given $C_{\text{ref}}$ by maximizing the likelihood:

[TABLE]

However, as the simulator is often very complex (e.g. biological simulators that try to mimic a 3D folding of the DNA given a set of constraints), the likelihood $p(C|\theta)$ may be intractable. As such, we directly target the posterior density $p(\theta|C_{\text{ref}})$ using data from the joint model ($\theta,C)\sim p(\theta)p(C|\theta)$, either via approximate Bayesian computation (Sequential Monte-Carlo ABC: SMC-ABC) or by estimating the posterior density with a conditional normalizing flow (Sequential neural posterior estimation: SNPE) snpea ; snpec .

SMC-ABC.

We use a variant of ABC coupled with sequential Monte-Carlo (SMC) smcabc . It consists of multiple rounds of ABC where, at each round, relevant $\{\theta^{k,*}\}_{k}$ are selected from the training set $\{(\theta^{n},C^{n})\}_{n}$ depending on a closeness criterion between $C$ and $C_{\text{ref}}$. We then associate weights $\{w^{k}\}_{k}$ to those selected $\{\theta^{k,*}\}_{k}$, and use the set $\{(\theta^{k,*},w^{k})\}_{k}$ to create the next population of $\{\theta^{n}\}_{n}$ for the next round of ABC. This sequential approach enables us to refine the relevant $\theta$ at each round. However, we need to define a metric for discriminating $(\theta^{n},\theta^{m})$ based on their associated observations $(C^{n},C^{m})$.

1. Metric 1** : Pearson correlation – ABC-Pearson.** To measure the closeness between $C$ and $C_{\text{ref}}$, the Pearson correlation is commonly used tjong ; nelly ; nelle . We find that the vector-based Pearson correlation averaged over all trans-contacts blocks is the most discriminative metric: each trans-contacts block of $C$ and $C_{\text{ref}}$ is vectorized and the Pearson correlation is computed between both. We then average all the correlations over the trans-contacts blocks (see Algorithm  2 in Appendix C.1). However, this metric is fine-tuned to this specific inference task.
2. Metric 2** : Data-driven summary statistic – ABC-CNN.** Instead of looking for a specific metric to compare $C$ to $C_{\text{ref}}$, we choose to use the classical $l^{2}$-norm. For this, we need a summary statistic $S$ that will extract the main features of $C$ and project it into a low-dimensional vector. One relevant candidate for a summary statistic is $\mathbb{E}\left[\theta|C\right]$ because with this one:

[TABLE]

where $\epsilon$ is the ABC-threshold.

When $\epsilon\to 0$, we don’t lose any first-order information when summarizing $C$ summstatabc . Moreover, $\mathbb{E}\left[\theta|C\right]$ is the analytical solution of the regression of $\theta$ on $C$ i.e.

[TABLE]

where $\mathcal{F}$ is the set of square integrable functions. As this statistic is unavailable, we learn it via a (deep) neural network (DNN) $S_{\phi}$ with parameters $\phi$ summstatabc . The DNN encoding $S_{\phi}$ is composed of a convolutional neural network (CNN) followed by a multi-layer perceptron (MLP). Using Monte Carlo estimator of (1) with $N$ samples $(\theta^{n},C^{n})\sim p(\theta)p(C|\theta)$, the DNN loss to be minimized in $\phi$ is then

[TABLE]

For large $N$, we expect $S_{\phi}(C)\approx\mathbb{E}\left[\theta|C\right]$. This approach has two phases: first learn $S_{\phi}$ by minimizing $\widehat{\mathcal{L}}_{\text{DNN}}(\phi)$, then run sequential ABC with this summary statistic and the $l^{2}$-norm as discriminative criterion (see Algorithm 3 in Appendix C.2).

SNPE – SBI-CNN.

SMC-ABC yields only samples from the target posterior distribution $p(\theta|C_{\text{ref}})$, but evaluating log-probabilities can be useful for downstream tasks. In contrast, a conditional normalizing flow $p_{\psi}(.|.)$ snpea ; snpec used to estimate the posterior distribution can both easily sample from the posterior and return the values of its log-probabilities. To ensure that $p_{\psi}(\theta|C_{\text{ref}})$ is close to $p(\theta|C_{\text{ref}})$, we minimize their Kullback–Leibler divergence ($D_{\mathrm{KL}}$) averaged over the observations $C$ as per

[TABLE]

After simplifications and using a Monte Carlo estimator, the flow is trained to minimize

[TABLE]

Once trained, we obtain an amortized estimator of the posterior densities $p(\theta|C)$ valid for any $C$. We just have to plug in $C_{\text{ref}}$ to get the estimated posterior density $p_{\psi}(.|C_{\text{ref}})$ (see Algorithm 4 in Appendix D). Since we are actually interested in the posterior at $C_{\text{ref}}$, parameters $\theta$ with very low posterior density may not be useful for learning $\psi$. Thus, we consider a sequential approach with several rounds of NPE to get an iterative refinement of the posterior estimate snpec . From the second round, $\theta^{n}$ are sampled from the latest estimated posterior found instead of the prior. This way, training samples are more informative about $C_{\text{ref}}$, gradually improving the learning of $\psi$. When the observations $C$ are high-dimensional (e.g. 2D-matrices), we encode them in a summary statistic $S$ using a convolutional neural network.

3 Numerical experiments

We showcase our methodology on two settings involving the genome of the yeast S. cerevisiae, for which we have access to the true position of all its centromeres: firstly, we run our inference pipeline on only S. cerevisiae’s first three chromosomes, secondly on its whole genome (16 chromosomes). We assess the performance of different inference methods by comparing their approximate posterior distributions to a ground-truth distribution consisting of Diracs on each dimension located at the true centromere positions. All experiments can be run on the CPU of a laptop, requiring $\sim 1$ h for the small genome and $\sim 5$ h for the entire genome.

3.1 Study case – small genome (3 chromosomes)

In low-dimensional settings, we can jointly infer $\theta$ given the entire contact map $C_{\text{ref}}$. For each inference method, we consider 11 rounds, each with a training dataset $\{(\theta^{n},C^{n})\}_{n}$ of size $10^{3}$. The summary statistic $S_{\phi}$ is pre-trained using a different training set of size $5\times 10^{3}$ with the optimizer Adam and a fixed learning rate $5\times 10^{-4}$. For SNPE, we use a masked autoregressive flow (MAF) and the version SNPE-C snpec from the Python package sbi sbipackage2024 (see Appendix E for more details).

3.2 High dimensional problem – whole genome (16 chromosomes)

When extending the analysis to the entire genome, we end up facing the curse of dimensionality: the space of parameters $\theta$ becomes too large to cover with few simulations and the contact maps $C$ are too big. As such, the resulting neural network encoding $S_{\phi}$ has too many parameters to be optimized. To reduce the dimension of the problem, we run $16$ parallel inferences (one per dimension of $\theta$) extracting each time only the informative part of the contact maps. With this approach, the space of $\theta$ is cut into several chromosome-length 1D intervals reducing the train set size. Let $C_{i}$ and $C_{\text{ref},i}$ be the $i^{\mathrm{th}}$ lines of blocks of matrices $C$ and $C_{\text{ref}}$, respectively. To infer $\theta_{i}$ with ABC-Pearson, we compute the vector-based Pearson correlation averaged over all blocks between $C_{i}$ and $C_{\text{ref},i}$. Concerning the data-driven summary statistic approaches ABC-CNN and SBI-CNN, the summary statistic $S_{\phi_{i}}$ also tries to project $C_{i}$ to $\theta_{i}$.

Using the redundancy of the data between rows of blocks, and to minimize the number of parameters of $\{S_{\phi_{i}}\}_{i}$, we consider a shared architecture where the CNN parameters are shared between chromosomes and the MLP ones are chromosome-specific. For SBI-CNN, we learn $16$ normalizing flows $p_{\psi_{i}}(\theta_{i}|S_{\phi_{i}}(C_{\text{ref},i}))$ (see Appendix F for more details and results).

4 Discussion

We present a novel methodology to infer the positions of the centromeres of the yeast S. cerevisiae using Hi-C contact maps. The probabilistic framework that we use allows us to quantify the uncertainty about the centromere candidates. Our entire inference pipeline is based on a large number of simulations relating centromere positions and contact maps. To mitigate computing bottlenecks, we have designed a simplified but efficient simulator that yields very convincing results when coupled with inferences on real experimental data.

In the case of a small genome, we obtained accurate inference of the centromere positions (Figure 1). The estimated densities for the summary statistic-based methods (ABC-CNN and SBI-CNN) are not very biased and peaky around the ground truth. In each dimension, $\theta$ is estimated at a precision under the resolution of the contact map $C_{\text{ref}}$ (Figure 6a) and the Euclidean distance to $\theta_{\text{ref}}$ is also under the resolution (Figure 1b). Moreover, SBI-CNN outperforms ABC-CNN that itself outperforms ABC-Pearson, reinforcing the use of a summary statistic and the flexibility of the normalizing flows. In the case of the whole genome, our approaches are not as accurate and could be improved. In some dimensions, $\theta$ is estimated at a precision under the resolution (Figure 7b), and we obtain peaky densities but in others the inference is not precise (Figure 2).

An advantage of our method is that we do not rely on any initialization or pre-localization: instead, we use an uninformative prior, setting each centromere randomly in the range of its corresponding chromosome. Also, our approach is naturally scalable: the pre-trained summary statistic could be reused for inference on centromeres of others yeasts without any re-training. To improve our approach, we will focus our efforts in developing a summary statistic independent of the size of the genome via notably the use of transformer architectures Vaswani2017 .

Acknowledgments and Disclosure of Funding

This work was supported by the ANR project BONSAI (grant ANR-23-CE23-0012-01) and by the ANR project Bayes-Duality (grant ANR-21-JSTM-0001).

Appendix A Contact maps

A contact map summarizes all the chromatin contacts observed over a population of DNA configurations. To construct it, we define the resolution of the map (the length of the chromosome fragment that will represent one pixel in the map). Each chromosome is then cut into fragments and each entry of the map represents the contact counts of any fragment with another over the population of DNA. This creates a matrix by blocks of interactions between chromosomes. Usually, we represent them by a heatmap.

Appendix B The simulator

The goal of the simulator is to create the upper trans-contact blocks of a contact map $C$ rapidly given the centromere positions $\theta$. We want to mimic the peak of interaction that appears in those blocks, as well as some rare interactions that can occur among the population of DNA.

Given the $L$ chromosome lengths in bp $\{l_{i}\}_{1\leq i\leq L}$, the centromere positions $\theta$ are sampled from the prior $\mathcal{U}(\underset{1\leq i\leq L}{\prod}[1,l_{i}-1])$. To create each contact map $C$, the process is described in Algorithm 1.

Our simulator outputs contacts maps that present some dissimilarities with Hi-C maps. If we compute the row-based averaged Pearson correlation between $C_{\text{ref}}$ and $C$ simulated from $\theta_{\text{ref}}$ as in [16], we get a correlation of $0.18$ in the case of $3$ chromosomes and $0.12$ in the case of $16$ chromosomes, which is quite low.

However, concerning the inference task, our simulator estimates the centromere positions $\theta$ nearly as well on synthetic data (Figure 5) as on Hi-C data (Figure 1).

Appendix C SMC-ABC

C.1 With the metric Pearson correlation – ABC-Pearson

One of the inference methods used is sequential ABC with the metric vector-based Pearson correlation averaged over all trans-contacts blocks.

When $\epsilon_{\text{corr}}\to 1$, $p(\theta|\text{corr}(C,C_{\text{ref}})\geq\epsilon_{\text{corr}})\to p(\theta|C_{\text{ref}})$.

C.2 With a summary statistic and the classical $l^{2}$-norm – ABC-CNN

The other ABC approach uses a pre-learned summary statistic $S_{\phi}$.

For $S_{\phi}$ informative enough, and when $\epsilon\to 0$,

[TABLE]

Appendix D SNPE – SBI-CNN

The last inference approach is SNPE based on normalizing flows and the pre-learned summary statistic $S_{\phi}$.

It is a sequential method: in the first round, $\theta$ is drawn from an uninformative prior. From the next rounds, $\theta$ is drawn from a proposal: the posterior density estimated from the previous round. This way, $\theta$ is more informative about $C_{\text{ref}}$ and the inference is expected to be refined across rounds.

Appendix E Small genome inference

We work with the S. cerevisiae’s first three chromosomes. $\theta$ is directly inferred from the entire contact map $C_{\text{ref}}$. We present a benchmark of metrics to assess the performance of the different inference methods : they evaluate both the proximity of the samples to $\theta_{\text{ref}}$ and the closeness of the densities to the ‘true’ posterior $\delta_{\theta_{\text{ref}}}$.

Appendix F Whole genome inference

To reduce the dimension of the problem, we carry 16 parallel inferences: one per dimension of $\theta$. Thus, we have 16 $1D$ inference problems where the parameter $\theta_{i}$ is drawn from a Uniform prior whose range is the size of the chromosome $i$ in bp. The simulator creates the $i^{\mathrm{th}}$ row of trans-contact blocks of a contact map $C$ (denoted $C_{i})$. All the inference methods target the posterior $p(\theta_{i}|C_{\text{ref},i})$. We need also to learn 16 summary statistics $\{S_{\phi_{i}}\}_{i}$ to project each row of trans-contact blocks $C_{i}$ to $\theta_{i}$.

$S_{\phi_{i}}$ is a CNN to capture the information of $C_{i}$ followed by an MLP to project this information into $\theta_{i}$. On the one hand, as the rows of trans-contact blocks $C_{i}$ are quite similar, we choose a shared architecture for the CNN between chromosomes. On the other hand, each MLP depends on the size of each chromosome so a chromosome-specific architecture is thus needed for this part of the DNN.

For the SBI method, we also need to learn $16$ normalizing flows. As for the 3-chromosomes case, we choose a MAF as well as SNPE-C for the experiments.