Simulation-based inference of yeast centromeres

Eloïse Touron

Univ. Grenoble Alpes, Inria
CNRS, Grenoble INP, LJK, France
eloise.touron@inria.fr

Pedro L. C. Rodrigues

Univ. Grenoble Alpes, Inria
CNRS, Grenoble INP, LJK, France
pedro.rodrigues@inria.fr

Julyan Arbel

Univ. Grenoble Alpes, Inria
CNRS, Grenoble INP, LJK, France
julyan.arbel@inria.fr

Nelle Varoquaux

TIMC, Univ. Grenoble Alpes
CNRS, Grenoble INP, France
nelle.varoquaux@univ-grenoble-alpes.fr

Michael Arbel

Univ. Grenoble Alpes, Inria
CNRS, Grenoble INP, LJK, France
michael.arbel@inria.fr

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.

# 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 [6, 20]. 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 investiga-
tion due to their essential role in many biological processes: they facilitate chromosome segregation
through the formation of the kinetochore [2] 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) [4] 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) [14] or chromatin
immunoprecipitation (ChIP) [12]. However, these approaches often remain imprecise and may even fail to infer centromeres for some species [7]. To bypass these limitations, newer methods have been proposed, which use Hi-C contact maps [13, 18]. 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_{\text{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. [5] to construct the reference Hi-C contact map. Centromere candidates $\theta = (\theta_1, \dots, \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 [1], 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 [10]. 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 [10]. 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-contact 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_{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_{ref}$ by maximizing the likelihood:

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

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_{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) [15, 8].

**SMC-ABC.** We use a variant of ABC coupled with sequential Monte-Carlo (SMC) [17]. 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_{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)$.

**Metric 1: Pearson correlation – ABC-Pearson.** To measure the closeness between $C$ and $C_{ref}$, the Pearson correlation is commonly used [16, 13, 18]. 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_{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.

**Metric 2: Data-driven summary statistic – ABC-CNN.** Instead of looking for a specific metric to compare $C$ to $C_{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}[\theta|C]$ because with this one:

$$ \mathbb{E}[\theta | \|S(C) - S(C_{ref})\| \le \epsilon] \xrightarrow{\epsilon \to 0} \mathbb{E}[\theta|C_{ref}], $$

where $\epsilon$ is the ABC-threshold.

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

$$ \mathbb{E}[\theta|C] = \underset{S \in \mathcal{F}}{\operatorname{argmin}} \mathbb{E}[\|S(C) - \theta\|_2^2], \quad (1) $$

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$ [11]. 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

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

For large $N$, we expect $S_{\phi}(C) \approx \mathbb{E}[\theta|C]$. This approach has two phases: first learn $S_{\phi}$ by minimizing $\hat{\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_{ref})$, but evaluating log-probabilities can be useful for downstream tasks. In contrast, a conditional normalizing flow $p_{\psi}(\cdot|\cdot)$ [15, 8] 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_{ref})$ is close to $p(\theta|C_{ref})$, we minimize their Kullback-Leibler divergence ($D_{KL}$) averaged over the observations $C$ as per

$$ \mathbb{E}_C[D_{KL}(p(\cdot|C)\|p_{\psi}(\cdot|C))]. $$
<figure><img src="image_1.png" /><figcaption>Figure 1: Inference using ABC-Pearson, ABC-CNN, and SBI-CNN (a). Color shades increase from lightest to darkest across rounds. Densities are estimated with the 5% best $\theta$ according to the ABC criterion or sampled from the flow. We also report the mean Euclidean distance between $\theta$ and $\theta_{ref}$, computed over the 5% best-performing samples in the top right corner (b). The horizontal dashed line stands for the resolution of the contact map $C_{ref}$ (in bp) in the top right figure. Results with SBI-CNN are uniformly better and both approaches based on data-driven summary statistics have errors smaller than the resolution of the contact maps.</figcaption></figure>

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

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

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_{ref}$ to get the estimated posterior density $p_{\psi}(\cdot|C_{ref})$ (see Algorithm 4 in Appendix D). Since we are actually interested in the posterior at $C_{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 [8]. 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_{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 ~ 1 h for the small genome and ~ 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_{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 [8] from the Python package sbi [3] (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
<figure><img src="image_2.png" /><figcaption>Figure 2: Inference using **ABC-Pearson**, **ABC-CNN**, and **SBI-CNN**. Color shades increase from lightest to darkest across rounds. Densities are estimated with the 5% best $\theta$ according to the ABC criterion or sampled from the flow. In some dimensions, the densities are very peaky and centered around $\theta_i$ (e.g. chromosome 4, 13, 15) but in others, the inference is not precise (e.g. chromosome 1, 6, 10). Data-driven summary statistics approaches do not outperform Pearson correlation-based method.</figcaption></figure>

$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_{ref,i}$ be the $i^{th}$ lines of blocks of matrices $C$ and $C_{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_{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_{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_{ref}$ (Figure 6a) and the Euclidean distance to $\theta_{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 [19].

## 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).

## References

[1] Jean-Michel Arbona, Sébastien Herbert, Emmanuelle Fabre, and Christophe Zimmer. Inference of the physical properties of yeast chromatin through Bayesian analysis of whole nucleus simulations. In *Genome Biology*, 2017.

[2] Kerry S. Bloom. Centromeric heterochromatin: the primordial segregation machine. *Annu. Rev. Genet.*, 48:457–484, 2014.

[3] Jan Boelts, Michael Deistler, Manuel Gloeckler, Álvaro Tejero-Cantero, Jan-Matthis Lueckmann, Guy Moss, Peter Steinbach, Thomas Moreau, Fabio Muratore, Julia Linhart, Conor Durkan, Julius Vetter, Benjamin Kurt Miller, Maternus Herold, Abolfazl Ziaeemehr, Matthijs Pals, Theo Gruner, Sebastian Bischoff, Nasty a Krouglova, Richard Gao, Janne K. Lappalainen, Bálint Mucsányi, Felix Pei, Auguste Schulz, Zinovia Stefanidi, Pedro Rodrigues, Cornelius Schröder, Faried Abu Zaid, Jonas Beck, Jaivardhan Kapoor, David S. Greenberg, Pedro J. Gonçalves, and Jakob H. Macke. sbi reloaded: a toolkit for simulation-based inference workflows, 2024.

[4] Guillaume Cottarel, James H. Shero, Philip Hieter, and Johannes H. Hegemann. A 125-base-pair CEN6 DNA fragment is sufficient for complete meiotic and mitotic centromere functions in Saccharomyces cerevisiae. *Mol Cell Biol*, 9(8):3342–3349, 1989.

[5] Zhijun Duan, Mirela Andronescu, Kevin Schutz, Sean McIlwain, Yoo Jung Kim, Choli Lee, Jay Shendure, Stanley Fields, Anthony Blau, and William S Noble. A three-dimensional model of the yeast genome. *Nature*, 465:363–367, 2010.

[6] Kyle P Eagen. Principles of chromosome architecture revealed by Hi-C. In *Trends Biochem Sci.*, 2018.

[7] Jonathan L. Gordon, Kevin P. Byme, and Kenneth H. Wolfe. Mechanisms of Chromosome Number Evolution in Yeast. *PLoS Genet*, 7, 2011.

[8] David S. Greenberg, Marcel Nonnenmacher, and Jakob H. Macke. Automatic posterior transformation for likelihood-free inference. In *ICML*, 2019.

[9] Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel two-sample test. In *Journal of Machine Learning Research*, 2012.

[10] Maxim Imakaev, Geoffrey Fudenberg, Rachel P. McCord, Natalia Naumova, Anton Goloborodko, Bryan Lajoie, Job Dekker, and Leonid Mirny. Iterative correction of Hi-C data reveals hallmarks of chromosome organization. *Nature Methods*, 9:999–1003, 2012.

[11] Bai Jiang, Tung yu Wu, Charles Zheng, and Wing H. Wong. Learning summary statistic for approximate Bayesian computation via deep neural network. In *Statistica Sinica*, 2018.
[12] Philippe Lefrançois, Ghia M. Euskirchen, Raymond K. Auerbach, Joel Rozowsky, Theodore Gibson, Christopher M. Yellman, Mark Gerstein, and Michael Snyder. Efficient yeast ChIP-Seq using multiplex short-read DNA sequencing. *BMC Genomics*, 10(37), 2009.

[13] Hervé Marie-Nelly, Martial Marbouty, Axel Cournac, Gianni Liti, Gilles Fischer, Christophe Zimmer, and Romain Koszul. Filling annotation gaps in yeast genomes using genome-wide contact maps. In *Bioinformatics*, 2014.

[14] Angela Nietzel, Mariano Rocchi, Heike Starke, Anita Heller, Wolfgang Fiedler, Iwona Wlodarska, Ivan Loncarevic, Volkmar Beensen, Uwe Claussen, and Thomas Liehr. A new multicolor-FISH approach for the characterization of marker chromosomes: centromere-specific multicolor-FISH (cenM-FISH). *Human Genetics*, 2001.

[15] George Papamakarios and Iain Murray. Fast e-free inference of simulation models with Bayesian conditional density estimation. In *NeurIPS*, 2016.

[16] Harianto Tjong, Ke Gong, Lin Chen, and Frank Alber. Physical tethering and volume exclusion determine higher-order genome organization in budding yeast. In *Genome Research*, 2012.

[17] Tina Toni, David Welch, Natalja Strelkowa, Andreas Ipsen, and Michael P.H Stumpf. Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. In *Journal of The Royal Society Interface*, 2008.

[18] Nelle Varoquaux, Ivan Liachko, Ferhat Ay, Joshua N Burton, Jay Shendure, Maitreya J Dunham, Jean-Philippe Vert, and William S Noble. Accurate identification of centromere locations in yeast genomes using Hi-C. In *Nucleic Acids Research*, 2015.

[19] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In *Advances in Neural Information Processing Systems*, 2017.

[20] Joachim Wolff, Rolf Backofen, and Björn Grüning. Loop detection using Hi-C data with HiCExplorer. In *GigaScience*, volume 11, 2022.
A Contact maps

A contact map summarizes all the chromatin contacts observed over a population of DNA configura-
tions. 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.

<figure><img src="image_3.png" /><figcaption>Figure 3: Process to construct a contact map in the case of 2 chromosomes.</figcaption></figure>

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 *θ*. 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≤*L*, the centromere positions *θ* are sampled from the prior *U*([1, *l*ᵢ - 1]). To create each contact map *C*, the process is described in Algorithm 1.

**Algorithm 1 Simulator of contact maps**

**Input:** *L* chromosome lengths in bp {lᵢ}₁≤i≤*L*, resolution of the contact map in bp *r* (e.g. *r* = 32 kb), centromere positions *θ*

**Return:** the upper trans-contact blocks of a simulated contact map *C* at the resolution *r* bp.

choose the size of the peaks of interaction: sample $\sigma^2$ from $U(0.1, 10)$

choose the intensity of interaction $\alpha$ to simulate the DNA population size: sample $\alpha$ from $\mathcal{U}(\llbracket 1, 1000 \rrbracket)$

for each chromosome pair (i, j) do

define a block of interaction $C_{ij}$ of size $(\frac{l_i}{r}, \frac{l_j}{r})$

define the center of the peak $(\theta_i, \theta_j)$

apply Gaussian density $N((\theta_i/r, \theta_j/r), \sigma^2)$ to the pixels of the block $C_{ij}$

multiply each pixel of $C_{ij}$ by the intensity factor $\alpha$

add Gaussian noise up to 10% of the maximal value of $C_{ij}$ to mimic the rare contacts:

construct a random matrix $M_{ij}$ of size $(\frac{l_i}{r}, \frac{l_j}{r})$ where each pixel is sampled from $N(\max(C_{ij}) \times 0.05, (\max(C_{ij}) \times 0.05)^2)$, then add $M_{ij}$ to $C_{ij}$

end for

return a simulated contact map *C* at resolution *r* bp
<figure><img src="image_4.png" /><figcaption>Figure 4: Hi-C map and our simulated map in the case of a small genome (resolution 32 kb).</figcaption></figure>

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).

<figure><img src="image_5.png" /><figcaption>Figure 5: Inference using **ABC-Pearson**, **ABC-CNN**, and **SBI-CNN** from synthetic data (a). Color shades increase from lightest to darkest across rounds. Densities are estimated with the 5% best $\theta$ according to the ABC criterion or sampled from the flow. We also report the mean Euclidean distance between $\theta$ and $\theta_{\text{ref}}$, computed over the 5% best-performing samples in the top right corner (b). The horizontal dashed line stands for the resolution of the contact map $C$ (in bp) in the top right figure. Results with data-driven summary statistics approaches are uniformly better even if all approaches have errors smaller than the resolution of the contact maps.</figcaption></figure>
# 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.

**Algorithm 2 SMC-ABC based on Pearson correlation inspired from [17]**

**Input:** T rounds, prior $\pi$, train set of size $N$, acceptance size $M$, perturbation kernel $K = \mathcal{N}(., \sigma^2\text{Id})$ ($\sigma = \text{resolution (bp)}$)

**Return:** $\theta \sim p(\theta|\text{corr}(C, C_{\text{ref}}) \ge \epsilon_{\text{corr}})$

**round** $t = 0$

- sample $\theta^n \sim \pi$, and $C^n \sim p(.,|\theta^n)$, $n \in [1, N]$

- compute $\text{corr}(C^n, C_{\text{ref}})$ and keep the top 5% of $\{\theta^n\}_n$ in terms of the highest correlation: $\{\theta^{m,0}, m \in [1, M]\}$

- compute weights $\{w^{m,0} = \frac{1}{M}, m \in [1, M]\}$

**output round** $t = 0$: $\{(\theta^{m,0}, w^{m,0})\}_{m \in [1, M]}$

**for** $0 < t < T$ **do**

**round** $t$

- from the previous accepted $\{\theta^{m,t-1}\}_{m \in [1, M]}$, sample $\{\bar{\theta}^k, k \in [1, M]\}$ from multinomial $M(\{\theta^{m,t-1}\}_m, \{w^{m,t-1}\}_m)$ with replacement

- perturb $\frac{N}{M}$ times the $M$ samples $\bar{\theta}^k$ to have $N$ samples $\theta^n$

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

- check that $\theta^n$ is in the prior bound otherwise, set $\theta^n \leftarrow \bar{\theta}^k$

- from this set $\{\theta^n\}_{n \in [1, N]}$, sample $C^n \sim p(.,|\theta^n)$, $n \in [1, N]$

- compute $\text{corr}(C^n, C_{\text{ref}})$ and keep the top 5% of $\{\theta^n\}_n$ in terms of the highest correlation: $\{\theta^{m,t}, m \in [1, M]\}$

- compute corresponding weights

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

**output round** $t$: $\{(\theta^{m,t}, w^{m,t})\}_{m \in [1, M]}$

**end for**

**return** accepted samples $\theta^n \sim p(\theta|\text{corr}(C^n, C_{\text{ref}}) \ge \epsilon_{\text{corr}})$

When $\epsilon_{\text{corr}} \to 1$, $p(\theta|\text{corr}(C, C_{\text{ref}}) \ge \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$.

**Algorithm 3** ABC with learned summary statistic inspired from [11]

**Input:** (deep) neural network (DNN) $S_\phi$, threshold $\epsilon$, Euclidean norm in $\mathbb{R}^n$, simulator, prior $p$

**Return:** Samples $\theta$ from the estimated posterior density $p(\cdot | \|S_\phi(C) - S_\phi(C_{\text{ref}})\| \le \epsilon)$

**Stage 1:** learn the summary statistic $S_\phi(\cdot)$ s.t. $S_\phi(C) \approx \mathbb{E}[\theta|C]$

generate a train set $(\theta^n, C^n)$ from $p(\theta)p(C|\theta)$

train a DNN $S_\phi$ on this train set with the loss to minimize in $\phi$

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

output $S_\phi(\cdot)$ s.t. $S_\phi(C) \approx \mathbb{E}[\theta|C]$

**Stage 2:** run ABC with the learned summary statistic $S_\phi$ and the criterion $\|S_\phi(C) - S_\phi(C_{\text{ref}})\| \le \epsilon$

**return** accepted samples $\theta^n \sim p(\cdot | \|S_\phi(C^n) - S_\phi(C_{\text{ref}})\| \le \epsilon)$

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

$$ p(\theta | \|S_\phi(C) - S_\phi(C_{\text{ref}})\| \le \epsilon) \to p(\theta|S_\phi(C_{\text{ref}})) \approx p(\theta|C_{\text{ref}}). $$

# 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.

**Algorithm 4** SNPE inspired from [15] and [8]

**Input:** T rounds, posterior density estimator $p_\psi$, simulator, prior $p$, simulation budget $N$, observation $C_{\text{ref}}$, pre-learned summary statistic $S_\phi$

**Return:** The estimated posterior density $p_\psi(\cdot|S_\phi(C_{\text{ref}}))$

**for** round $t = 1, ..., T$ **do**

**if** $t = 1$ **then** $p_t = p$

**end if**

**for** $n = 1, ..., N$ **do**

sample $\theta^n \sim p_t$

sample $C^n \sim p(\cdot|\theta^n)$

**end for**

train the posterior estimator $p_\psi$ on $\mathcal{D} = \{(\theta^n, C^n)\}_n$ with the loss to minimize in $\psi$

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

use $p_\psi$ to construct the estimated posterior: $p_\psi(\cdot|S_\phi(C_{\text{ref}}))$.

define the proposal for the next round: $p_t(\theta) = p_\psi(\theta|S_\phi(C_{\text{ref}}))$

**end for**

**return** samples $\theta^n \sim p_\psi(\theta|S_\phi(C_{\text{ref}}))$
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}}}$.

<figure><img src="image_6.png" /><figcaption>Figure 6: We report the absolute error per dimension of $\theta$ between the mean computed over the 5% best-performing samples and $\theta_{\text{ref}}$ (a) as well as the Maximum Mean Discrepancy (MMD) (b) and the Wasserstein-2 distance (c) between $p(\theta|C_{\text{ref}})$ and $\delta_{\theta_{\text{ref}}}$ [9].</figcaption></figure>

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^{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.
<figure><img src="image_7.png" /><figcaption>Figure 7: We report the mean Euclidean distance between $\theta$ and $\theta_{\text{ref}}$ (a), computed over the 5% best-performing samples, the absolute error per dimension of $\theta$ between the mean $\theta$ computed over the 5% best-performing samples and $\theta_{\text{ref}}$ (b) as well as the MMD (c) and the Wasserstein-2 distance (d) between $p(\theta|C_{\text{ref}})$ and $\delta_{\theta_{\text{ref}}}$. The horizontal dotted line stands for the resolution of the contact map $C_{\text{ref}}$ (in bp) in the top figures.</figcaption></figure>