remarkRemark hypothesisHypothesis claimClaim 3D Genome Organization in Diploid
OrganismsA. Belyaeva, K. Kubjas, L. J. Sun and C. Uhler

# Identifying 3D Genome Organization in Diploid Organisms via Euclidean
Distance Geometry††thanks: Submitted. Anastasiya Belyaeva was supported by an
NSF Graduate Research Fellowship (1122374), the Abdul Latif Jameel World Water
and Food Security Lab (J-WAFS) at MIT and the MIT J-Clinic for Machine
Learning and Health. Kaie Kubjas was supported by the European Union’s Horizon
2020 research and innovation programme: Marie Skłodowska-Curie grant agreement
No. 748354, research carried out at LIDS, MIT and Team PolSys, LIP6, Sorbonne
University. Caroline Uhler was partially supported by NSF (DMS-1651995), ONR
(N00014-17-1-2147 and N00014-18-1-2765), IBM, and a Simons Investigator Award.

Anastasiya Belyaeva Laboratory for Information and Decision Systems,
Department of Electrical Engineering and Computer Science, and Institute for
Data, Systems and Society, Massachusetts Institute of Technology, Cambridge,
MA<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Kaie Kubjas Department
of Mathematics and Systems Analysis, Aalto University, Espoo, Finland ().
<EMAIL_ADDRESS>Lawrence J. Sun22footnotemark: 2 Caroline
Uhler22footnotemark: 2

###### Abstract

The spatial organization of the DNA in the cell nucleus plays an important
role for gene regulation, DNA replication, and genomic integrity. Through the
development of chromosome conformation capture experiments (such as 3C, 4C,
Hi-C) it is now possible to obtain the contact frequencies of the DNA at the
whole-genome level. In this paper, we study the problem of reconstructing the
3D organization of the genome from such whole-genome contact frequencies. A
standard approach is to transform the contact frequencies into noisy distance
measurements and then apply semidefinite programming (SDP) formulations to
obtain the 3D configuration. However, neglected in such reconstructions is the
fact that most eukaryotes including humans are diploid and therefore contain
two copies of each genomic locus. We prove that the 3D organization of the DNA
is not identifiable from distance measurements derived from contact
frequencies in diploid organisms. In fact, there are infinitely many solutions
even in the noise-free setting. We then discuss various additional
biologically relevant and experimentally measurable constraints (including
distances between neighboring genomic loci and higher-order interactions) and
prove identifiability under these conditions. Furthermore, we provide SDP
formulations for computing the 3D embedding of the DNA with these additional
constraints and show that we can recover the true 3D embedding with high
accuracy from both noiseless and noisy measurements. Finally, we apply our
algorithm to real pairwise and higher-order contact frequency data and show
that we can recover known genome organization patterns.

###### keywords:

3D genome organization; Hi-C; diploid organisms; Euclidean distance geometry;
semidefinite programming, systems of polynomial equations.

51K05, 92E10, 90C22, 52C25, 14P05.

##  1 Introduction

It is now well established that the spatial organization of the genome in the
cell nucleus plays an important role for cellular processes including gene
regulation, DNA replication, and the maintenance of genomic integrity [11, 39,
40]. Notably, a recent study [43] showed a causal link between three-
dimensional (3D) genome organization and gene regulation, where gene
repositioning was induced and subsequent changes in gene expression were
observed. This motivates the development of methods to reconstruct the 3D
structure of the genome to study its functions.

The genetic information in cells is contained in the DNA, which is organized
into chromosomes and packed into the cell nucleus. Chromosome confirmation
capture techniques (such as 3C, 4C, Hi-C, Capture-C) have enabled the
interrogation of the contact frequencies between pairs of genomic loci at the
whole-genome scale [12, 37, 24, 20]. In Hi-C, for example, interacting
chromosome regions are crosslinked (i.e., frozen), the DNA is then fragmented,
the crosslinked fragments are ligated, and paired-end sequencing is applied to
the ligation products and mapped to a reference genome [24]. By binning the
genome and ascribing each read pair into the corresponding bin, one obtains a
contact frequency matrix between genomic loci that is commonly of the size
$10^{6}\times 10^{6}$.

Different computational approaches for reconstructing the 3D genome
organization from contact frequency data have been considered. Distance-based
approaches convert contact frequencies $F_{ij}$ into spatial distances
$D_{ij}$ and find a Euclidean embedding of the points in 3D [14, 46, 23, 35].
Ensemble methods such as MCMC5C and BACH [36, 19] learn a set of possible 3D
structures by defining a probabilistic model for contact frequencies and
generating an ensemble of structures via MCMC sampling. Other ensemble methods
include molecular dynamics simulations that model DNA as a polymer and output
an ensemble of 3D structures [25, 27, 13, 32]. Finally, statistical methods
have also been proposed that directly model contact counts instead of
distances, using for example the Poisson distribution [42], and maximize the
log-likelihood of the data to infer the 3D genome organization.

Almost all existing methods make the simplifying assumption that the genome is
haploid, when in fact most organisms of interest including humans are diploid,
i.e. there are two copies of each chromosome known as _homologous
chromosomes_. For example, human cells contain two copies of 23 chromosomes
each. The challenge is that the contact frequency data from chromosome
conformation capture experiments is generally _unphased_ , meaning that the
copies of each chromosome cannot be distinguished. As a result, if the DNA is
modeled as a string of beads containing two copies of each bead $i$ for $1\leq
i\leq n$, then the measured contact frequencies result in an $n\times n$
matrix, from which we would like to infer the 3D embedding of $2n$ points.
This problem cannot be solved by classical methods for 3D genome
reconstruction methods such as those mentioned above. With significant
experimental efforts, phased data can be obtained [7, 8] and used in order to
reconstruct the 3D genome organization [4]. However, such data is rare and
costly.

Figure 1: Schematic of the diploid genome. Nucleus with green, blue and red
curves depicting three homologous pairs of chromosomes. In the unphased
setting, the measured distance between loci $i$ and $j$ corresponds to the sum
of the four distances (denoted in purple) between two pairs of homologous loci
$x_{i},y_{i}$ and $x_{j},y_{j}$.

In this paper, we provide a computational method for inferring the 3D diploid
organization of the genome without relying on phased data. In particular, we
consider a distance-based approach and use Euclidean distance geometry to
obtain the 3D diploid structure of the genome. The precise mathematical
problem considered in this paper is as follows and illustrated in Figure 1.
DNA is modeled as a string of beads, that contains two copies of each bead $i$
for $1\leq i\leq n$. We would like to infer the location of the two copies of
each bead, which we denote by $x_{i}\in\mathbb{R}^{3}$ and
$y_{i}\in\mathbb{R}^{3}$. Since for unphased data, the two copies of each bead
cannot be distinguished, the problem is to identify the 3D configuration
($2n\times 3$ matrix), i.e. $x_{1},\dots x_{n},y_{1},\dots
y_{n}\in\mathbb{R}^{3}$ (up to translation and rotation), from the composite
distance measurements $D_{ij}$, $1\leq i\neq j\leq n$ ($n\times n$ matrix),
corresponding to the sum of the distances between either copy of bead $i$ and
$j$, i.e.,

$D_{ij}=\|x_{i}-x_{j}\|^{2}+\|x_{i}-y_{j}\|^{2}+\|y_{i}-x_{j}\|^{2}+\|y_{i}-y_{j}\|^{2}.$

In the haploid or phased setting, this problem boils down to the standard
Euclidean distance geometry problem. This problem has a long history: in the
classical setting with no missing values, this problem can be solved via the
classical multidimensional scaling (cMDS) algorithm that is based on spectral
decomposition followed by dimensionality reduction; see [9] for an overview.
Other approaches for the Euclidean embedding and completion problems,
including in the presence of missing values, are non-convex formulations [15,
28] as well as semidefinite relaxations [1, 16, 5, 26, 44, 45].

A naive approach in the unphased diploid setting is to assume that the four
distances that make up our measured composite distance $D_{ij}$ are equal and
solve the corresponding Euclidean embedding problem. However, it is evident
from single-cell imaging studies that the four distances in $D_{ij}$ can be
wildly different [3, 30]. Hence this approach cannot provide realistic
embeddings. While, a simple dimension argument ($6n$ variables versus
$\binom{n}{2}$ constraints) suggests that the 3D genome configuration is
uniquely identifiable, one of the main results of our paper is that the 3D
diploid genome configuration is not identifiable from unphased data. In fact,
we show that there are infinitely many configurations that satisfy the
constraints imposed by $D_{ij}$, even in the noiseless setting (Section 2,
Theorem 2.1).

We therefore consider additional biologically relevant and experimentally
measurable constraints and study identifiability of the 3D diploid structure
under these constraints. First, we take into account distances between
neighboring beads, i.e. $\|x_{i}-x_{i+1}\|^{2}$ and $\|y_{i}-y_{i+1}\|^{2}$ on
each chromosome. While we show that this yields unique identifiability for
configurations in 2D, there are still infinitely many configurations in 3D,
which is of primary interest for genome modeling (Section 3, Proposition 3.1,
Proposition 3.3). To obtain identifiability in 3D, we consider adding
constraints based on contact frequencies between three or more loci
simultaneously. The measurement of such higher-order contact frequencies has
recently been enabled by experimental assays such as SPRITE [33], C-walks [31]
and GAM [2]. We prove that this information can be used to uniquely identify
the 3D genome organization from unphased data in the noiseless setting
(Section 4, Theorem 4.1).

Finally, we provide an SDP formulation for obtaining the 3D diploid
configuration from noisy measurements (Section 5) and show based on simulated
data that our algorithm has good performance and that it is able to recover
known genome organization patterns when applied to real contact frequency data
collected from human lymphoblastoid cells (Section 6).

##  2 Unidentifiability from pairwise distance constraints

In the remainder of the paper we denote the true but unknown coordinates of
the homologous loci by $x^{*}_{i}$ and $y^{*}_{i}$ and the corresponding
noiseless distances by $D^{*}_{ij}$ while the symbols $x_{i}$ and $y_{i}$
denote the variables that we want to solve for. While from a biological
perspective the relevant setting is when $x_{i},y_{i}\in\mathbb{R}^{3}$,
results that hold more generally will be stated in $\mathbb{R}^{d}$. The main
result of this section is Theorem 2.1, which characterizes the set of
solutions given by the constraints $D^{*}_{ij}$ in dimension $d\leq 3$. In
particular, it establishes non-identifiability of the 3D genome structure from
pairwise distance measurements in the diploid unphased setting.

######  Theorem 2.1.

Let $d\leq 3$ and $n\geq 2d+3$. Then
$(x_{1},\ldots,x_{n},y_{1},\ldots,y_{n})\in(\mathbb{R}^{d})^{2n}$ satisfies

(1)
$D^{*}_{ij}=\|x_{i}-x_{j}\|^{2}+\|x_{i}-y_{j}\|^{2}+\|y_{i}-x_{j}\|^{2}+\|y_{i}-y_{j}\|^{2}\;\text{
for all }\,1\leq i\neq j\leq n$

if and only if it satisfies

(2) $x_{i}+y_{i}=x^{*}_{i}+y^{*}_{i}\;\text{ and
}\;\|x_{i}\|^{2}+\|y_{i}\|^{2}=\|x^{*}_{i}\|^{2}+\|y^{*}_{i}\|^{2}\;\text{ for
all }\,1\leq i\leq n$

up to translations and rotations in $\mathbb{R}^{d}$ and permutations of
$x_{i}$ and $y_{i}$.

As a consequence, the measurements $D_{ij}^{*}$ identify the location of each
pair of homologous loci $(x_{i},y_{i})$ up to a sphere with center
$(x^{*}_{i}+y^{*}_{i})/2$ and radius $\|x^{*}_{i}-y^{*}_{i}\|/2$. Namely, the
points $x_{i},y_{i}$ lie opposite to each other anywhere on this sphere.
Unless $x^{*}_{i}=y^{*}_{i}$ for all $i$, i.e., all spheres have radius $0$,
this set is infinite in dimensions $d>1$ and hence the configuration is
unidentifiable.

In the remainder of this section, we will prove Theorem 2.1. The two
inclusions in Theorem 2.1 are proven in Lemma 2.2 and Lemma 2.6. In Lemma 2.4
it is shown that the distance $\|x_{i}-y_{i}\|$ within each homologous pair is
fixed given the pairwise distances $D^{*}_{ij}$. This result is used to prove
Lemma 2.6.

######  Lemma 2.2.

Let $(x_{1},\ldots,x_{n},y_{1},\ldots,y_{n})\in(\mathbb{R}^{d})^{2n}$ satisfy

(3) $x_{i}+y_{i}=x^{*}_{i}+y^{*}_{i}\text{ and
}\|x_{i}\|^{2}+\|y_{i}\|^{2}=\|x^{*}_{i}\|^{2}+\|y^{*}_{i}\|^{2}\;\text{ for
all }\,1\leq i\leq n.$

Then

$\|x_{i}-x_{j}\|^{2}+\|x_{i}-y_{j}\|^{2}+\|y_{i}-x_{j}\|^{2}+\|y_{i}-y_{j}\|^{2}=D^{*}_{ij}\;\text{
for all }\,1\leq i\neq j\leq n.$

######  Proof 2.3.

Observe that for each pair $x_{i},y_{i}$ satisfying the equations Eq. 3, it
holds that

$\displaystyle D^{*}_{ij}$ $\displaystyle=$ $\displaystyle
2\cdot(\|x^{*}_{i}\|^{2}+\|y^{*}_{i}\|^{2})+2\cdot(\|x^{*}_{j}\|^{2}+\|y^{*}_{j}\|^{2})-2(x^{*}_{i}+y^{*}_{i})\cdot(x^{*}_{j}+y^{*}_{j})$
$\displaystyle=$ $\displaystyle
2\cdot(\|x_{i}\|^{2}+\|y_{i}\|^{2})+2\cdot(\|x_{j}\|^{2}+\|y_{j}\|^{2})-2(x_{i}+y_{i})\cdot(x_{j}+y_{j})$
$\displaystyle=$
$\displaystyle\|x_{i}-x_{j}\|^{2}+\|x_{i}-y_{j}\|^{2}+\|y_{i}-x_{j}\|^{2}+\|y_{i}-y_{j}\|^{2}.$

This completes the proof.

Next we will show that the distance between homologous pairs is uniquely
determined by the $D^{*}_{ij}$

######  Lemma 2.4.

Let $d\leq 3$ and $n\geq 2d+3$. Then for each $1\leq i\leq n$ the quantity
$\|x_{i}-y_{i}\|$ is identifiable from the constraints imposed by the
$D^{*}_{ij}$, i.e., for any solution
$(x_{1},\ldots,x_{n},y_{1},\ldots,y_{n})\in(\mathbb{R}^{d})^{2n}$ to the
equations defined by the $D^{*}_{ij}$ in Eq. 1, the quantity $\|x_{i}-y_{i}\|$
is constant.

The constraint $d\leq 3$ is due to our proof technique. The condition $n\geq
2d+3$ is necessary for unique identifiability of the distance between
homologous pairs of loci.

######  Proof 2.5.

Without loss of generality we assume that $i=1$ and show that
$\|x_{1}-y_{1}\|$ is equal to some constant. First, we perform a shift on the
solution so that $x_{1}=-y_{1}=v$. Since shifts preserve distances, they in
particular preserve the equality constraints Eq. 1. Hence,

$D^{*}_{1j}=\|v-x_{j}\|^{2}+\|v-y_{j}\|^{2}+\|-v-x_{j}\|^{2}+\|-v-y_{j}\|^{2}.$

Expanding this out into dot products and simplifying yields

$D^{*}_{1j}=4\|v\|^{2}+2(\|x_{j}\|^{2}+\|y_{j}\|^{2}).$

Let $j\neq k$ be both not equal to $1$. Then substituting the above leads to

$D^{*}_{1j}+D^{*}_{1k}-D^{*}_{jk}=8\|v\|^{2}+2(x_{j}+y_{j})\cdot(x_{k}+y_{k}).$

Defining $T_{jk}:=D^{*}_{1j}+D^{*}_{1k}-D^{*}_{jk}$ and
$s_{j}:=\sqrt{2}(x_{j}+y_{j})$, this is equivalent to

$T_{jk}-8\|v\|^{2}=s_{j}\cdot s_{k}.$

Let $T^{\prime}$ be the $(d+1)\times(d+1)$ submatrix of $T$ satisfying
$T^{\prime}_{ij}=T_{i+1,j+d+2}$, i.e. the rows of $T^{\prime}$ correspond to
the rows $2,3,\ldots,d+2$ of $T$ and the columns of $T^{\prime}$ correspond to
the columns $d+3,d+4,\ldots,2d+3$ of $T$. We now show that for generic
configurations $\det(T^{\prime})\neq 0$. Since $\det(T^{\prime})$ can be
written as a polynomial in the coordinates $x_{i}$ and $y_{i}$, then
$\det(T^{\prime})\neq 0$ for generic configurations as long as it does not
identically vanish. Hence it suffices to present one configuration where
$\det(T^{\prime})$ is nonzero. For $d\leq 3$ we can check this using random
configurations.

Since $T^{\prime}$ has full rank, then the matrix determinant lemma implies
that

(4)
$\det(T^{\prime}-8J\|v\|^{2})=(1-8\|v\|^{2}\mathds{1}^{T}(T^{\prime})^{-1}\mathds{1})\det(T^{\prime}),$

where $\mathds{1}$ denotes the all ones vector. Note that the scalar
$\mathds{1}^{T}T^{\prime-1}\mathds{1}$ is fixed and $(\det T^{\prime})\neq 0$.
Furthermore, since $T^{\prime}-8J\|v\|^{2}$ is formed from the dot products
between $d$-dimensional vectors, it has rank at most $d$ and therefore
$\det(T^{\prime}-8J\|v\|^{2})=0$ due to $T^{\prime}-8J\|v\|^{2}$ being a
$(d+1)\times(d+1)$ matrix. Hence,
$(1-8\|v\|^{2}\mathds{1}^{T}(T^{\prime})^{-1}\mathds{1})\det(T^{\prime})=0$,
which is a linear equation in terms of $\|v\|^{2}$. As a consequence, it has a
unique solution for $\|v\|^{2}$ and thus the distance between the homologous
pair $x_{1},y_{1}$ is fixed as long as $n\geq 2d+3$.

We next characterize all solutions to the constraints imposed by the
$D^{*}_{ij}$.

######  Lemma 2.6.

Let $d\leq 3$ and $n\geq 2d+3$. Let
$(x_{1},\ldots,x_{n},y_{1},\ldots,y_{n})\in(\mathbb{R}^{d})^{2n}$ be a
solution to

$\|x_{i}-x_{j}\|^{2}+\|x_{i}-y_{j}\|^{2}+\|y_{i}-x_{j}\|^{2}+\|y_{i}-y_{j}\|^{2}=D^{*}_{ij}\;\text{
for all }\,1\leq i\neq j\leq n.$

Then

$x_{i}+y_{i}=x^{*}_{i}+y^{*}_{i}\text{ and
}\|x_{i}\|^{2}+\|y_{i}\|^{2}=\|x^{*}_{i}\|^{2}+\|y^{*}_{i}\|^{2}\;\text{ for
all }\,1\leq i\leq n$

up to translations and rotations in $\mathbb{R}^{d}$ and permutations of
$x_{i}$ and $y_{i}$.

######  Proof 2.7.

Without loss of generality we perform a translation on the solution such that
$x_{1}=-y_{1}=v$ for some vector $v$. By Lemma 2.4 the quantity
$\|x_{k}-y_{k}\|$ is constant for each $1\leq k\leq n$ and thus also $\|v\|$
is constant. Since for any $j\neq 1$ it holds that
$D^{*}_{1j}=4\|v\|^{2}+2(\|x_{j}\|^{2}+\|y_{j}\|^{2})$, also
$\|x_{j}\|^{2}+\|y_{j}\|^{2}$ is constant and hence
$\|x_{i}\|^{2}+\|y_{i}\|^{2}=\|x^{*}_{i}\|^{2}+\|y^{*}_{i}\|^{2}$ for all
$1\leq i\leq n$.

Similarly to the proof of Lemma 2.4, if we define
$T_{jk}=D^{*}_{1j}+D^{*}_{1k}-D^{*}_{jk}$ and $s_{j}=\sqrt{2}(x_{j}+y_{j})$,
we find that

$T_{jk}-8\|v\|^{2}=s_{j}\cdot s_{k}.$

Because we have access to the diagonal constraints now, this relationship
holds for all $j,k$ and not just $j\neq k$. Thus $T-8J\|v\|^{2}$ is a
symmetric $(n-1)\times(n-1)$ matrix admitting a rank $d$ factorization. Let
$S$ be the matrix formed with the vectors $s_{j}$. We then have
$T-8J\|v\|^{2}=SS^{T}$. There is a result on rank factorizations of symmetric
matrices that any other factorization $T-8J\|v\|^{2}=S^{\prime}S^{\prime T}$
satisfies $S=S^{\prime}Q$ for some orthogonal matrix $Q$ [22, Proposition
3.2]. Thus for any other solution $s^{\prime}_{j}$, we have
$s_{j}=s^{\prime}_{j}Q$, implying all solutions are simply orthogonal
transformations of each other (rotations, reflections, etc.)

In summary, we have shown that once we have fixed $x_{1}+y_{1}=0$ via
translation, then the quantities $x_{j}+y_{j}$ are unique up to orthogonal
transformations and the quantities $\|x_{j}\|^{2}+\|y_{j}\|^{2}$ are unique.

##  3 Distance constraints between neighboring loci

In Section 2, we showed that the 3D genome configuration is not identifiable
from pairwise distance constraints available from typical (unphased) contact
frequency maps. In order to gain identifiability, we next consider adding
other biological constraints to the problem formulation that are generally
available or can be measured. In particular, since DNA can be viewed as a
string of connected beads, we use the distance between adjacent beads as an
additional constraint. The distance between neighboring beads can be derived
empirically for example from imaging studies [29, 21]; see also our
experimental results in Section 6. The additional mathematical constraints
are:

$\|x_{i}-x_{i+1}\|=\|x^{*}_{i}-x^{*}_{i+1}\|\text{ and
}\|y_{i}-y_{i+1}\|=\|y^{*}_{i}-y^{*}_{i+1}\|\text{ for }1\leq i\leq n-1,$

where $x^{*}_{1},x^{*}_{2},\ldots,x^{*}_{n}$ and
$y^{*}_{1},y^{*}_{2},\ldots,y^{*}_{n}$ correspond to consecutive beads on
homologous chromosomes; see Figure 2.

In this section we show the following results: under the additional distance
constraints between neighboring loci, we prove that identifiability can be
obtained in the 2D setting (Proposition 3.1). However, in the 3D setting we
prove that there are still infinitely many 3D configurations even with these
additional distance constraints (Proposition 3.3).

For the proofs of Proposition 3.1 and Proposition 3.3 we recall from Theorem
2.1 that $(x_{i},y_{i})$ and $(x^{*}_{i},y^{*}_{i})$ are diametrically
opposite points on the same sphere. Denote the $i$-th sphere by $S_{i}$ and
let it have center $c_{i}$ and radius $r_{i}$. Then $\|c_{i}-x_{i}\|=r_{i}$
and $2c_{i}-x_{i}=y_{i}$.

Figure 2: Distance constraints between neighboring beads. Green and blue
curves depict two homologous pairs of chromosomes. For the green curves
distances between neighboring genomic regions are shown by black lines.

######  Proposition 3.1.

For $n\geq 3$ and generic
$(x^{*}_{1},\ldots,x^{*}_{n},y^{*}_{1},\ldots,y^{*}_{n})\in(\mathbb{R}^{2})^{2n}$,
there is a unique point
$(x_{1},\ldots,x_{n},y_{1},\ldots,y_{n})\in(\mathbb{R}^{2})^{2n}$ satisfying
the equations

(5) $\begin{gathered}x_{i}+y_{i}=x^{*}_{i}+y^{*}_{i}\text{ and
}\|x_{i}\|^{2}+\|y_{i}\|^{2}=\|x^{*}_{i}\|^{2}+\|y^{*}_{i}\|^{2}\text{ for
}1\leq i\leq n,\\\ \|x_{i}-x_{i+1}\|=\|x^{*}_{i}-x^{*}_{i+1}\|\text{ and
}\|y_{i}-y_{i+1}\|=\|y^{*}_{i}-y^{*}_{i+1}\|\text{ for }1\leq i\leq
n-1.\end{gathered}$

######  Proof 3.2.

We have $y_{1}=2c_{1}-x_{1}$ and $y_{2}=2c_{2}-x_{2}$. Plugging this into
$\|y_{1}-y_{2}\|=\|y^{*}_{1}-y^{*}_{2}\|$ gives

$\displaystyle\|y^{*}_{1}-y^{*}_{2}\|$ $\displaystyle=$
$\displaystyle\|(2c_{1}-x_{1})-(2c_{2}-x_{2})\|^{2}$ $\displaystyle=$
$\displaystyle\|(2c_{1}-2c_{2})-(x_{1}-x_{2})\|^{2}$ $\displaystyle=$
$\displaystyle\|2c_{1}-2c_{2}\|^{2}+\|x_{1}-x_{2}\|^{2}-2(2c_{1}-2c_{2})\cdot(x_{1}-x_{2}).$

The quantities $\|2c_{1}-2c_{2}\|^{2}$ and $\|x_{1}-x_{2}\|^{2}$ are fixed.
This implies that the quantity $(2c_{1}-2c_{2})\cdot(x_{1}-x_{2})$ is fixed.
Since we know $\|x_{1}-x_{2}\|$ and $c_{1}\neq c_{2}$ holds by genericness,
then there are two possible angles for $x_{1}-x_{2}$ (this is where we use the
2D constraint) and thus that there are two possible solutions for
$x_{1}-x_{2}$.

Because $x_{1},x_{2}$ are constrained to lie on circles, the solutions for
$x_{1}$ are the intersection points of the first circle and the second circle
translated by $x_{1}-x_{2}$ and the solutions for $x_{2}$ are the intersection
points of the second circle and the first circle translated by $x_{2}-x_{1}$.
Hence each solution for $x_{1}-x_{2}$ leads to at most two possible solutions
for $(x_{1},x_{2})$. In turn this implies there are at most four solutions for
$x_{2}$.

Figure 3: Identifiability in the 2D setting with neighboring distance
constraints. Two solutions for $x_{2}$ are obtained by translating the circle
centered at $c_{1}$ by $x^{*}_{2}-x^{*}_{1}$ (this new circle is colored blue)
and intersecting it with the circle centered at $c_{2}$. The other two
solutions are obtained by reflecting the blue circle over the line through
$c_{1}$ and $c_{2}$ (this new circle is colored green) and intersecting it
with the circle centered at $c_{2}$. The true solution for $x_{2}$ is colored
black and the three alternative solutions for $x_{2}$ are colored red.

We now investigate the four solutions. The first two solutions are obtained by
translating the circle centered at $c_{1}$ by $x^{*}_{2}-x^{*}_{1}$ and
intersecting it with the circle centered at $c_{2}$, see Figure 3. One of the
two solutions is $x^{*}_{2}$. The other two solutions are reflections of these
two solutions over the line from $c_{1}$ to $c_{2}$.

Let $x^{*}_{1},x^{*}_{2},c_{1},c_{2}$ be fixed. They determine four possible
solutions for $x_{2}$. We will show that these four solutions are different
from the four solutions we get from considering
$x^{*}_{2},x^{*}_{3},c_{2},c_{3}$ for generic $x^{*}_{3},c_{3}$ (apart from
$x^{*}_{2}$).

If either of the reflected solutions over the line from $c_{2}$ to $c_{3}$
coincides with one of the four original solutions, then we can perturb $c_{3}$
away from the line from $c_{2}$ to $c_{3}$ to change these solutions. If the
solution that is the intersection point of the circle centered at $c_{2}$ and
the translation by $x^{*}_{2}-x^{*}_{3}$ of the circle centered at $c_{3}$
(different from $x^{*}_{2}$) coincides with one of the four original
solutions, then we can perturb $x^{*}_{3}$. This changes $x^{*}_{2}-x^{*}_{3}$
and hence the second intersection point of the circle centered at $c_{2}$ and
the translation by $x^{*}_{2}-x^{*}_{3}$ of the circle centered at $c_{3}$.

A similar argument can be used to show that $x_{3},\ldots,x_{n-1}$ have unique
solutions. Given a unique solution for $x_{2}$, there are two solutions for
$x_{1}$ if and only if $x^{*}_{2}$ lies on the line from $c_{1}$ to $c_{2}$.
This is however not a generic configuration. A similar argument applies for
$x_{n}$.

Despite having uniqueness in 2D, we do not have uniqueness in 3D as shown in
the following proposition.

######  Proposition 3.3.

Let $n\in\mathbb{N}$. For generic
$(x^{*}_{1},\ldots,x^{*}_{n},y^{*}_{1},\ldots,y^{*}_{n})\in(\mathbb{R}^{3})^{2n}$,
there are infinitely many points
$(x_{1},\ldots,x_{n},y_{1},\ldots,y_{n})\in(\mathbb{R}^{3})^{2n}$ satisfying
equations Eq. 5.

######  Proof 3.4.

If $n=1$, then $x_{1}^{*},y_{1}^{*}$ can be chosen randomly with the
constraint that $x_{1}^{*}\neq y_{1}^{*}$. Then $x_{1}$ and $y_{1}$ can be any
points on the sphere $S_{1}$ defined by $x_{1}^{*},y_{1}^{*}$. Now assume that
$n\geq 2$. Fix any $x_{1}^{*},y_{1}^{*}$ such that $x_{1}^{*}\neq y_{1}^{*}$.
Choose two circles $C_{1}$ and $C^{\prime}_{1}$ on the sphere $S_{1}$ defined
by $x_{1}^{*},y_{1}^{*}$ that intersect at two points one of which is
$x_{1}^{*}$. The circle $C_{1}$ is the intersection of $S_{1}$ and another
sphere $T_{1}$. Let $x_{2}^{*}$ be the center of the sphere $T_{1}$. Let
$C^{\prime\prime}_{1}$ be the circle on $S_{1}$ that consists of points
antipodal to $C^{\prime}_{1}$. Then $C^{\prime\prime}_{1}$ is also an
intersection of $S_{1}$ and another sphere $T^{\prime\prime}_{2}$. Let
$y_{2}^{*}$ be the center of the sphere $T^{\prime\prime}_{2}$. We use the
same procedure to construct $x_{3}^{*}$ and $y_{3}^{*}$ from $x_{2}^{*}$ and
$y_{2}^{*}$, $x_{4}^{*}$ and $y_{4}^{*}$ from $x_{3}^{*}$ and $y_{3}^{*}$ etc.

The only condition on $x_{1}^{*}$ and $y_{1}^{*}$ is $x_{1}^{*}\neq
y_{1}^{*}$, hence $(x_{1}^{*},y_{1}^{*})$ is a generic point in
$\mathbb{R}^{3}\times\mathbb{R}^{3}$. The condition that $C_{1}$ is a circle
on the sphere $S_{1}$ containing $x_{1}^{*}$ is equivalent to $x_{2}^{*}$
being any point in $\mathbb{R}^{3}$ outside the line through $x_{1}^{*}$ and
$y_{1}^{*}$. Similarly, the condition that $C^{\prime}_{1}$ is a circle on the
sphere $S_{1}$ containing $x_{1}^{*}$ is equivalent to $y_{2}^{*}$ being any
point in $\mathbb{R}^{3}$ outside the line through $x_{1}^{*}$ and
$y_{1}^{*}$. The condition that $C_{1}$ and $C^{\prime}_{1}$ intersect at two
different points of $S_{1}$ is equivalent to the normal vector of the tangent
plane of $S_{1}$ at $x_{1}^{*}$ and the normal vectors of the planes defined
by $C_{1}$ and $C^{\prime}_{1}$ being linearly independent. Hence
$(x_{1}^{*},x_{2}^{*},y_{1}^{*},y_{2}^{*})$ is a generic point in
$(\mathbb{R}^{3})^{4}$. Similar arguments can be used to show that
$(x^{*}_{1},x^{*}_{2},\ldots,x^{*}_{n},y^{*}_{1},y^{*}_{2},\ldots,y^{*}_{n})$
is a generic point in $(\mathbb{R}^{3})^{2n}$.

Now consider points $x_{n}$ and $y_{n}$ in an $\varepsilon$-neighborhood of
$x_{n}^{*}$ and $y_{n}^{*}$. Consider the spheres that are centered at $x_{n}$
and $y_{n}$ and have radii $\|x^{*}_{n-1}-x^{*}_{n}\|\text{ and
}\|y^{*}_{n-1}-y^{*}_{n}\|$. The intersections of these spheres with $S_{n-1}$
give circles $\tilde{C}_{n-1}$ and $\tilde{C}^{\prime\prime}_{n-1}$ that are
perturbations of circles $C_{n-1}$ and $C^{\prime\prime}_{n-1}$. In
particular, the intersection of the circle $\tilde{C}_{n-1}$ and the circle
$\tilde{C}^{\prime}_{n-1}$ that consists of points antipodal to
$\tilde{C}^{\prime\prime}_{n-1}$ consists of two points for $\varepsilon$
small enough. Choosing $x_{n-1}$ to be the intersection point corresponding to
$x_{n-1}^{*}$ and $y_{n-1}$ its antipodal gives points $x_{n-1},y_{n-1}$
satisfying $\|x_{n-1}-x_{n}\|=\|x^{*}_{n-1}-x^{*}_{n}\|\text{ and
}\|y_{n-1}-y_{n}\|=\|y^{*}_{n-1}-y^{*}_{n}\|$.

Assuming that $\varepsilon$ is small enough, then $x_{n-1}$ and $y_{n-1}$ are
in small neighborhoods of $x_{n-1}^{*}$ and $y_{n-1}^{*}$, and we can continue
the same procedure to find $x_{n-2}$ and $y_{n-2}$ from $x_{n-1}$ and
$y_{n-1}$, $x_{n-3}$ and $y_{n-3}$ from $x_{n-2}$ and $y_{n-2}$ etc. In
particular, we can find $x_{1},\ldots,x_{n-1},y_{1},\ldots,y_{n-1}$ satisfying
equations Eq. 5 for every $x_{n}$ and $y_{n}$ in an $\varepsilon$-neighborhood
of $x_{n}^{*}$ and $y_{n}^{*}$.

The previous proposition suggests that there are two degrees of freedom for
choosing $x_{1},\ldots,x_{n},y_{1},\ldots,y_{n}$ on each homologous pair and
thus that finite identifiability requires two additional algebraically
independent constraints per homologous pair. Similarly this suggests that
unique identifiability requires three additional algebraically independent
constraints per homologous pair, where each endpoint of a chromosome needs to
be included in at least one of the additional constraints.

##  4 Identifiability from higher-order contact constraints

In Section 3, we showed that considering distances between neighboring beads
only yields identifiability in 2D but not in 3D. In the following, we consider
adding further constraints that are becoming widely available from
experimental data, namely higher-order contact frequencies between three or
more loci as measured by experimental assays such as SPRITE [33], C-walks [31]
and GAM [2]. We express these constraints mathematically by letting
$F\in\mathbb{R}^{m\times m\times\cdots\times m}$ be a contact frequency
tensor, where $F_{x_{i_{1}},x_{i_{2}},\ldots,x_{i_{k}}}$ measures the contact
frequency between loci $i_{1},i_{2},\ldots,i_{k}$ with coordinates
$x_{i_{1}},x_{i_{2}},\ldots,x_{i_{k}}$. In the unphased setting, we can only
measure a combination of contact frequencies over the homologous loci
$\\{x_{i_{1}},y_{i_{1}}\\}\times\\{x_{i_{2}},y_{i_{2}}\\}\times\ldots\times\\{x_{i_{k}},y_{i_{k}}\\}$,
which we denote by $F_{i_{1}i_{2}\ldots i_{k}}$. In addition, as for 2-way
interactions, we turn contact frequencies into “distances” by defining
$D_{i_{1}i_{2}\ldots i_{k}}:=1/F_{i_{1}i_{2}\ldots i_{k}}$.

(a)

(b)

Figure 4: Higher-order distance constraints. (a) Three loci
$x_{i_{1}},x_{i_{2}},x_{i_{3}}$, located on the same chromosome, are depicted.
In the phased setting, the higher-order distance
$D_{x_{i_{1}}x_{i_{2}}x_{i_{3}}}$ is defined as the sum of the distances (pink
dashed lines) of the three loci $x_{i_{1}},x_{i_{2}},x_{i_{3}}$ to their
centroid (pink point). Green and blue curves depict two different chromosomes.
(b) Illustrates the definition of $D_{i_{1}i_{2}i_{3}}$ in the unphased
setting. Green, blue and red curves depict neighborhoods around three
homologous loci $(x_{i_{1}},y_{i_{1}}),(x_{i_{2}},y_{i_{2}})$ and
$(x_{i_{3}},y_{i_{3}})$. From these homologous loci 8 possible higher-order
distances can be defined (colored dashed lines) based on the 8 centroids
depicted in the figure. The higher-order distance $D_{i_{1}i_{2}i_{3}}$ is
defined as the minimum of these 8 distances (achieved here by the three black
dashed line segments).

In the following, we provide our interpretation of distances in the higher-
order setting. For simplicity, we first describe the higher-order distances
for three loci in the phased setting. Since $F_{x_{i_{1}}x_{i_{2}}x_{i_{3}}}$
counts how often the three loci come together, we interpret
$D_{x_{i_{1}}x_{i_{2}}x_{i_{3}}}$ as the sum of the distances of the three
loci $x_{i_{1}},x_{i_{2}},x_{i_{3}}$ to their centroid (Figure 4a). We next
provide a generalization to the unphased setting. For three homologous loci
$(x_{i_{1}},y_{i_{1}}),(x_{i_{2}},y_{i_{2}})$ and $(x_{i_{3}},y_{i_{3}})$,
their contact frequency can be formed by 8 possible triples, namely
$(x_{i_{1}},x_{i_{2}},x_{i_{3}})$, $(x_{i_{1}},x_{i_{2}},y_{i_{3}})$,
$(x_{i_{1}},y_{i_{2}},x_{i_{3}})$, $(y_{i_{1}},x_{i_{2}},x_{i_{3}})$,
$(x_{i_{1}},y_{i_{2}},y_{i_{3}})$, $(y_{i_{1}},y_{i_{2}},x_{i_{3}})$,
$(y_{i_{1}},x_{i_{2}},y_{i_{3}})$, and $(y_{i_{1}},y_{i_{2}},y_{i_{3}})$. We
will assume that one of the triples constitutes the majority of the observed
contact frequency count and hence we define $D_{i_{1}i_{2}i_{3}}$ as the
minimum over all 8 higher order distances. This is illustrated in Figure 4b.
Generalizing from three to $k$ loci, our higher-order distance definition then
becomes

$D_{i_{1}i_{2}\ldots
i_{k}}=\min_{z_{i_{j}}\in\\{x_{i_{j}},y_{i_{j}}\\}}\left(\sum_{j=1}^{k}\|z_{i_{j}}-(z_{i_{1}}+\ldots+z_{i_{k}})/k\|^{2}\right).$

In the following, we prove our main result; namely we show that the distance
constraints of order 3 (3-way distances) together with the previously
considered pairwise distance constraints and distance constraints among
consecutive beads results in unique identifiability of the 3D genome
configuration (Theorem 4.1). In fact, only very few order 3 distance
constraints are required for unique identifiability. As we show in Theorem 4.1
it is sufficient that the first and last bead of each chromosome be contained
in an order 3 distance constraint. This is a reasonable constraint given that
methods such as SPRITE, C-walks and GAM measure higher-order interactions over
the whole genome. These insights are of interest experimentally since they
suggest that the methods can restrict the measurement of such higher order
constraints to first and last beads of each chromosome, known as telomeres.

######  Theorem 4.1.

Let $m$ be the number of chromosome pairs, let $n_{1},n_{2},\ldots,n_{m}$ be
the number of domains on chromosomes $1,2,\ldots,m$ and define
$n=n_{1}+n_{2}+\ldots+n_{m}$. Let $I\subseteq[n]\times[n]\times[n]$ be such
that each of
$1,n_{1},n_{1}+1,n_{1}+n_{2},\ldots,n_{1}+n_{2}+\ldots+n_{m-1}+1,n$ (labels of
domains at the beginning and at the end of each chromosome) is contained in at
least one triple in $I$. Let
$x^{*}_{1},\ldots,x^{*}_{n},y^{*}_{1},\ldots,y^{*}_{n}\in\mathbb{R}^{3}$ be
fixed such that

$\begin{gathered}\min_{z^{*}_{i}\in\\{x^{*}_{i},y^{*}_{i}\\}\text{ for
}i=k_{1},k_{2},k_{3}}\left(\sum_{j\in\\{k_{1},k_{2},k_{3}\\}}\|z^{*}_{j}-(z^{*}_{k_{1}}+z^{*}_{k_{2}}+z^{*}_{k_{3}})/3\|^{2}\right)=0\\\
\text{ for }(k_{1},k_{2},k_{3})\in I.\end{gathered}$

Consider the polynomial system:

(6) $\begin{gathered}x_{i}+y_{i}=x^{*}_{i}+y^{*}_{i}\text{ and
}\|x_{i}\|^{2}+\|y_{i}\|^{2}=\|x^{*}_{i}\|^{2}+\|y^{*}_{i}\|^{2}\text{ for
}1\leq i\leq n,\\\ \|x_{i}-x_{i+1}\|=\|x^{*}_{i}-x^{*}_{i+1}\|\text{ and
}\|y_{i}-y_{i+1}\|=\|y^{*}_{i}-y^{*}_{i+1}\|\\\ \text{ for
}i\in[n]\backslash\\{n_{1},n_{1}+n_{2},\ldots,n\\},\\\
\min_{z_{i}\in\\{x_{i},y_{i}\\}\text{ for
}i=k_{1},k_{2},k_{3}}\left(\sum_{j\in\\{k_{1},k_{2},k_{3}\\}}\|z_{j}-(z_{k_{1}}+z_{k_{2}}+z_{k_{3}})/3\|^{2}\right)=0\\\
\text{ for }(k_{1},k_{2},k_{3})\in I.\end{gathered}$

Then for generic $x^{*}_{1},\ldots,x^{*}_{n},y^{*}_{1},\ldots,y^{*}_{n}$, this
system has a unique solution in $(\mathbb{R}^{3})^{2n}$.

To prove Theorem 4.1, we will need two lemmas. Lemma 4.2 states that for a
fixed solution $(x_{1}^{*},y_{1}^{*})$ on a sphere $S_{1}$ and given distances
between solutions on $S_{1}$ and $S_{2}$, there are finitely many solutions
$(x_{2},y_{2})$ on the sphere $S_{2}$. Lemma 4.4 is an extension of Lemma 4.2.
It states that if one has finitely many solutions on a sphere $S_{i}$, then
given distances between neighboring beads, there are finitely many solutions
on any sphere connected to $S_{i}$.

######  Lemma 4.2.

Let $x^{*}_{1},x^{*}_{2},y^{*}_{1},y^{*}_{2}\in\mathbb{R}^{3}$ be fixed.
Consider the polynomial system:

(7)
$\begin{gathered}x_{2}+y_{2}=x^{*}_{2}+y^{*}_{2},\|x_{2}\|^{2}+\|y_{2}\|^{2}=\|x^{*}_{2}\|^{2}+\|y^{*}_{2}\|^{2},\\\
\|x^{*}_{1}-x_{2}\|=\|x^{*}_{1}-x^{*}_{2}\|\text{ and
}\|y^{*}_{1}-y_{2}\|=\|y^{*}_{1}-y^{*}_{2}\|.\end{gathered}$

For generic $x^{*}_{1},x^{*}_{2},y^{*}_{1},y^{*}_{2}$, this system has
finitely many solutions in $(\mathbb{R}^{3})^{2n}$.

######  Proof 4.3.

The first two equations of Eq. 7 say that $x_{2},y_{2}$ and
$x_{2}^{*},y_{2}^{*}$ are pairs of antipodal points on the same sphere. We
denote this sphere by $S_{2}$. The third equation says that $x_{2}$ is the
same distance from $x_{1}^{*}$ as $x_{2}^{*}$ is from $x_{1}^{*}$. Hence
$x_{2}$ must lie on the circle $C_{x_{2}}$ that is the intersection of $S_{2}$
and the sphere centered at $x_{1}^{*}$ and with radius
$\|x^{*}_{1}-x^{*}_{2}\|$. The last equation says that $y_{2}$ must lie on the
circle $C_{y_{2}}$ that is the intersection of $S_{2}$ and the sphere centered
at $y_{1}^{*}$ with radius $\|y^{*}_{1}-y^{*}_{2}\|$. We consider the circle
$C^{\prime}_{x_{2}}$ that consists of antipodal points to the circle
$C_{y_{2}}$ on the sphere $S_{2}$. The intersection of the circles $C_{x_{2}}$
and $C^{\prime}_{x_{2}}$ gives the solutions for $x_{2}$. Unless the two
circles are equal, they intersect at at most two points. Since $y_{2}$ is
antipodal to $x_{2}$, then for each $x_{2}$ there is a unique $y_{2}$. The
circles coincide if and only if $x_{1}^{*},y_{1}^{*}$ and the center of
$S_{2}$ are collinear.

######  Lemma 4.4.

Let $x^{*}_{1},\ldots,x^{*}_{n},y^{*}_{1},\ldots,y^{*}_{n}\in\mathbb{R}^{3}$
be fixed. Consider the polynomial system:

$\begin{gathered}x_{i}+y_{i}=x^{*}_{i}+y^{*}_{i}\text{ and
}\|x_{i}\|^{2}+\|y_{i}\|^{2}=\|x^{*}_{i}\|^{2}+\|y^{*}_{i}\|^{2}\text{ for
}2\leq i\leq n,\\\
\|x^{*}_{1}-x_{2}\|=\|x^{*}_{1}-x^{*}_{2}\|,\|y^{*}_{1}-y_{2}\|=\|y^{*}_{1}-y^{*}_{2}\|,\\\
\|x_{i}-x_{i+1}\|=\|x^{*}_{i}-x^{*}_{i+1}\|\text{ and
}\|y_{i}-y_{i+1}\|=\|y^{*}_{i}-y^{*}_{i+1}\|\text{ for }2\leq i\leq
n-1.\end{gathered}$

For generic $x^{*}_{1},\ldots,x^{*}_{n},y^{*}_{1},\ldots,y^{*}_{n}$, this
system has finitely many solutions in $(\mathbb{R}^{3})^{2n-2}$.

######  Proof 4.5.

By Lemma 4.2, there are finitely many antipodal pairs
$(x_{2},y_{2})\in\mathbb{R}^{3}\times\mathbb{R}^{3}$ on $S_{2}$ such that
$\|x^{*}_{1}-x_{2}\|=\|x^{*}_{1}-x^{*}_{2}\|$ and
$\|y^{*}_{1}-y_{2}\|=\|y^{*}_{1}-y^{*}_{2}\|$. Similarly, for each of these
antipodal pairs $(x_{2},y_{2})\in\mathbb{R}^{3}\times\mathbb{R}^{3}$ on
$S_{2}$, there are finitely many antipodal pairs
$(x_{3},y_{3})\in\mathbb{R}^{3}\times\mathbb{R}^{3}$ on $S_{3}$ satisfying
$\|x_{2}-x_{3}\|=\|x^{*}_{2}-x^{*}_{3}\|$ and
$\|y_{2}-y_{3}\|=\|y^{*}_{2}-y^{*}_{3}\|$ etc.

######  Proof 4.6 (Proof of Theorem 4.1).

We recall that the first line of the polynomial system Eq. 6 gives that
$x_{i},y_{i}$ are antipodal points on a sphere $S_{i}$. Consider a triple
$(k_{1},k_{2},k_{3})\in I$ that contains $1$ and the equation on the last line
of the polynomial system Eq. 6 corresponding to this triple. This equation
gives that $z_{k_{1}},z_{k_{2}},z_{k_{3}}$, where $z_{i}\in\\{x_{i},y_{i}\\}$,
coincide. Hence $z_{k_{1}},z_{k_{2}},z_{k_{3}}$ lie on the intersection of
$S_{k_{1}},S_{k_{2}},S_{k_{3}}$. Generically, if the intersection of three
spheres is non-empty in $\mathbb{R}^{3}$, then it consists of two points $P$
and $P^{\prime}$. This gives four possible solutions for $x_{1},y_{1}$: the
points $P,P^{\prime}$ and their antipodals on $S_{1}$. By Lemma 4.4, there are
finitely many solutions for $x_{2},\ldots,x_{n_{1}},y_{2},\ldots,y_{n_{1}}$
given these fixed solutions $x_{1},y_{1}$ on $S_{1}$. In the next two
paragraphs we will show that generically these finitely many solutions do not
contain antipodal points on any of the spheres $S_{2},\ldots,S_{n_{1}}$.

If there are two antipodal solutions on $S_{i}$, then we may assume that they
come either from the same solution on $S_{1}$ or antipodal solutions on
$S_{1}$, because we can perturb $S_{k_{1}},S_{k_{2}},S_{k_{3}}$ slightly to
change the other pair of solutions. First we will show that generically a
solution for $x_{i}$ on $S_{i}$ does not give a pair of antipodal solutions
for $x_{i+1}$ on $S_{i+1}$. If this was the case, then both the solution for
$x_{i}$ and its antipodal would have to lie on the plane that is perpendicular
to the line through the antipodal pair of solutions for $x_{i+1}$ on
$S_{i+1}$. This plane contains the centers of $S_{i}$ and $S_{i+1}$. Hence for
a solution for $x_{i}$, there is only one antipodal pair on solutions on
$S_{i+1}$. Thus for a generic distance between the solutions on $S_{i}$ and
$S_{i+1}$, a solution on $S_{i}$ does not give an antipodal pair of solutions
on $S_{n+1}$.

Secondly, suppose that two different solutions on $S_{i}$ give a pair of
antipodal solutions on $S_{i+1}$. We will show that when we perturb the
distance between solutions on $S_{i}$ and $S_{i+1}$, then we do not get an
antipodal pair anymore. Let $x_{i}$ and $x^{\prime}_{i}$ be two different
solutions on $S_{i}$ that give solutions $x_{i+1}$ and $2c_{i+1}-x_{i+1}$ on
$S_{i+1}$. Hence
$\|2c_{i+1}-x_{i+1}-x^{\prime}_{i}\|^{2}=\|x_{i+1}-x_{i}\|^{2}$. We want to
show that generically

$\|2c_{i+1}-(x_{i+1}+\epsilon)-x^{\prime}_{i}\|^{2}\neq\|x_{i+1}+\epsilon-
x_{i}\|^{2},$

where $x_{i+1}+\epsilon$ is the perturbed solution. Indeed, using the identity
$\|x_{i+1}-x_{i}\|^{2}=\|2c_{i+1}-x_{i+1}-x^{\prime}_{i}\|^{2}$ gives

$\|2c_{i+1}-(x_{i+1}+\epsilon)-x^{\prime}_{i}\|^{2}-\|x_{i+1}+\epsilon-
x_{i}\|^{2}=2\epsilon(x_{i}+x^{\prime}_{i}-2c_{i+1}).$

This quantity is equal to zero if and only if $\epsilon=0$ or $c_{i+1}$ is the
middle point of the line segment from $x_{i}$ to $x^{\prime}_{i}$. This is
generically not the case.

Using a triple $(k^{\prime}_{1},k^{\prime}_{2},k^{\prime}_{3})\in I$
containing $n_{1}$ and the equation for this triple, we get four possible
solutions for $x_{n_{1}},y_{n_{1}}$. Generically, only one of them coincides
with the finitely many solutions on $S_{n_{1}}$ that we get from the solutions
on $S_{1}$, because perturbing the spheres slightly (with keeping the
coinciding points fixed) perturbs the second intersection point of the three
spheres and we know that generically the finitely many points do not contain
antipodal points.

The unique solution on $S_{n}$ comes from one solution on each of the spheres
$S_{1},\ldots,S_{n_{1}-1}$: If this was not the case then two different
solutions on $S_{i}$ give the same solution on $S_{i+1}$. By the proof of
Proposition 3.1, the dot product $(c_{i}-c_{i+1})\cdot(x_{i}-x_{i+1})$ is
fixed. Hence for a fixed $x_{i+1}$, all possible solutions for $x_{i}$ lie on
a hyperplane and this hyperplane is perpendicular to $c_{i}-c_{i+1}$.
Therefore, if two solutions on $S_{i}$ give the same solution on $S_{i+1}$,
then they lie on a hyperplane perpendicular to $c_{i}-c_{i+1}$. By slightly
perturbing the sphere $S_{i+1}$, this is not the case anymore, and hence
generically a solution on $S_{i+1}$ comes from a unique solution on $S_{i}$.

##  5 Algorithms and implementation

So far, we derived a theoretical framework to establish when we have unique
and finite identifiability of the 3D configuration in the noiseless setting.
However, a unique solution does not necessarily mean that we can find it
efficiently, as in many cases finding the solution may be NP-hard. In
addition, we have so far not yet considered the noisy setting. In this
section, we show how to construct an optimization formulation to determine the
3D configuration efficiently.

We frame the 3D reconstruction problem as a Euclidean embedding problem, where
the coordinates $x_{1},\dots x_{n},y_{1},\dots y_{n}\in\mathbb{R}^{3}$ are
inferred from distances. Similar to ChromSDE [46], we formulate all distances
in terms of entries in the Gram matrix $G$, which tracks the dot products
between the $2n$ genomic regions. Namely, letting the column/row $i$ of $G$
correspond to $x_{i}$ and the column/row $n+i$ correspond to its homologous
locus $y_{i}$, then the distances are given by
$\|x_{i}-x_{j}\|^{2}=G_{i,i}+G_{j,j}-2G_{i,j}$,
$\|x_{i}-y_{j}\|^{2}=G_{i,i}+G_{n+j,n+j}-2G_{i,n+j}$ and
$\|y_{i}-y_{j}\|^{2}=G_{n+i,n+i}+G_{n+j,n+j}-2G_{n+i,n+j}$. It is natural to
work with the Gram matrix $G$, since it is rotation invariant. By imposing the
constraint $\sum_{i,j}G_{i,j}=0$ we can also fix the translational axis. Also
the additional distance constraints that we introduced in the previous
sections (LABEL:{theorem:non_identifiability}, Proposition 3.1, Theorem 4.1)
can be represented as linear constraints in terms of entries in $G$ as
follows:

  * •

Pairwise distances:

$g_{ij}(G):=G_{i,i}+G_{j,j}+G_{n+i,n+i}+G_{n+j,n+j}-G_{i,j}-G_{n+i,j}-G_{i,n+j}-G_{n+i,n+j}$

  * •

Distances between homologous pairs:

$g_{ii}(G):=G_{i,i}+G_{n+i,n+i}-2G_{i,n+i}$

  * •

Distances between neighboring beads:

$g_{i+}(G):=G_{i,i}+G_{i+1,i+1}-2G_{i,i+1}$

  * •

Distances of order 3 (can be generalized to higher orders):

$g_{ijk}(G):=\min_{l}(g_{ijkl}:l=1,\ldots,8)$

where

$\displaystyle g_{ijk1}(G):=$ $\displaystyle
G_{i,i}+G_{j,j}+G_{k,k}-G_{i,j}-G_{i,k}-G_{j,k},$ $\displaystyle
g_{ijk2}(G):=$ $\displaystyle
G_{i,i}+G_{j,j}+G_{n+k,n+k}-G_{i,j}-G_{i,n+k}-G_{j,n+k},$ $\displaystyle
g_{ijk3}(G):=$ $\displaystyle
G_{i,i}+G_{n+j,n+j}+G_{k,k}-G_{i,n+j}-G_{i,k}-G_{n+j,k},$ $\displaystyle
g_{ijk4}(G):=$ $\displaystyle
G_{i,i}+G_{n+j,n+j}+G_{n+k,n+k}-G_{i,n+j}-G_{i,n+k}-G_{n+j,n+k},$
$\displaystyle g_{ijk5}(G):=$ $\displaystyle
G_{n+i,n+i}+G_{j,j}+G_{k,k}-G_{n+i,j}-G_{n+i,k}-G_{j,k},$ $\displaystyle
g_{ijk6}(G):=$ $\displaystyle
G_{n+i,n+i}+G_{j,j}+G_{n+k,n+k}-G_{n+i,j}-G_{n+i,n+k}-G_{j,n+k},$
$\displaystyle g_{ijk7}(G):=$ $\displaystyle
G_{n+i,n+i}+G_{n+j,n+j}+G_{k,k}-G_{n+i,n+j}-G_{n+i,k}-G_{n+j,k},$
$\displaystyle g_{ijk8}(G):=$ $\displaystyle
G_{n+i,n+i}+G_{n+j,n+j}+G_{n+k,n+k}$ $\displaystyle-
G_{n+i,n+j}-G_{n+i,n+k}-G_{n+j,n+k}.$

Our objective is to determine a rank 3 solution of $G$, satisfying the above
constraints. However, this optimization problem is non-convex due to the rank
constraint, and we instead consider the standard relaxation: we minimize the
trace of the Gram matrix as an approximation to matrix rank [17]. The
resulting optimization problem then becomes the following semidefinite program
(SDP):

(8) $\displaystyle\underset{G}{\text{minimize}}$
$\displaystyle\operatorname{tr}(G)$ subject to $\displaystyle
g_{ii}(G)=D^{*}_{ii},\;1\leq i\leq n,$ $\displaystyle
g_{ij}(G)=D^{*}_{ij},\;1\leq i<j\leq n,$ $\displaystyle
g_{i+}(G)=D^{*}_{i+},\;i\in\Omega_{1},$ $\displaystyle
g_{ijk}(G)=D^{*}_{ijk},\;(i,j,k)\in\Omega_{2},$ $\displaystyle\sum_{1\leq
i,j\leq 2n}G_{i,j}=0,$ $\displaystyle G\succeq 0.$

Here, $D^{*}_{ii}$ denote the distances between homologous pairs computed from
the pairwise distances using Lemma 2.4, $D^{*}_{ij}$ denote the pairwise
distances, $D^{*}_{i+}$ denote the distances between neighboring beads, and
$D^{*}_{ijk}$ denote the distances between three loci (while one could also
consider 4 or higher order distance constraints, in our implementation we only
used 3-way distance constraints since higher-order contacts are extremely
sparse). The index set
$\Omega_{1}=[2n]\backslash\\{n_{1},n_{1}+n_{2},\ldots,n,n+n_{1},n+n_{1}+n_{2},\ldots,2n\\}$
corresponds to all beads that are not the last bead on a chromosome. The index
set $\Omega_{2}\subseteq[n]^{3}$ corresponds to all triples of beads with non-
zero contact frequencies.

In the noisy setting, which is relevant for biological data, we replace the
equality constraints by penalties in the loss function. Namely, using $D^{*}$
for the noiseless and $D$ for the noisy distances, we replace the equality
constraints of the form $g(G)=D^{*}$ by adding $(g(G)-D)^{2}$ to the objective
function. For the higher-order distance constraints of the form
$D^{*}_{ijk}=\min(g_{ijk1}(G),\ldots,g_{ijk8}(G))$ for $(i,j,k)\in\Omega_{2}$
we use slack variables and a convex relaxation using an atomic norm that
combines the $\ell_{2}$\- and $\ell_{1}$-norms. More precisely, we propose the
use of the following transformation in the noisy setting,

$D_{ijk}+\lambda_{ijkl}=g_{ijkl}(G)+s_{ijkl}\text{ for }l=1,2,\ldots,8,$

where $\lambda_{ijkl},s_{ijkl}\geq 0$ for all $i,j,k,l$ act as slack
variables. In general, for each triple $(i,j,k)$ we want one of the
$\lambda_{ijkl}$ to be close to $0$ and the sum over all $s_{ijkl}$ to be
small. Naively this can be done by placing $\sum s_{ijkl}+\sum\lambda_{ijkl}$
into the objective function. However, this would not enforce for each
$(i,j,k)$ at least one $\lambda_{ijkl}$ to be close to $0$. Instead we propose
to use

$\sum_{(i,j,k)\in\Omega_{2},1\leq l\leq
8}s_{ijkl}+\sqrt{\sum_{(i,j,k)\in\Omega_{2}}\left(\sum_{1\leq l\leq
8}\lambda_{ijkl}\right)^{2}}.$

The $\ell_{2}$-norm will push down the $\sum_{l}\lambda_{ijkl}$ for each
$(i,j,k)$, while the $\ell_{1}$ norm will drive at least one of the
$\lambda_{ijkl}$ to zero, which is precisely the desired behavior. The
quantity $\sqrt{\sum_{i,j,k}\left(\sum_{l}\lambda_{ijkl}\right)^{2}}$ is an
atomic norm as defined in [6] with the set of atoms

$\begin{gathered}\mathcal{A}=\\{(\lambda_{ijkl}):\sum_{i,j,k}\left(\sum_{l}\lambda_{ijkl}\right)^{2}=1\text{
and }\\\
\qquad\qquad\qquad\qquad\qquad\sum_{i,j,k}\lambda_{ijkl_{ijk}}^{2}=1\text{ for
}l_{ijk}=1,\ldots,8,(i,j,k)\in\Omega_{2}\\}.\end{gathered}$

Then the optimization problem in the noisy setting becomes:

(9) $\displaystyle\underset{G,s,\lambda}{\text{minimize}}$
$\displaystyle\rho\operatorname{tr}(G)+\sum_{1\leq i\leq
n}(g_{ii}(G)-D_{ii})^{2}+\sum_{1\leq i<j\leq n}(g_{ij}(G)-D_{ij})^{2}$
$\displaystyle+\sum_{i\in\Omega_{1}}(g_{i+}(G)-D_{i+})^{2}+\sum_{(i,j,k)\in\Omega_{2},1\leq
l\leq 8}s_{ijkl}$
$\displaystyle+\sqrt{\sum_{(i,j,k)\in\Omega_{2}}\left(\sum_{1\leq l\leq
8}\lambda_{ijkl}\right)^{2}}$ subject to $\displaystyle
D_{ijk}+\lambda_{ijkl}=g_{ijkl}(G)+s_{ijkl},\;(i,j,k)\in\Omega_{2},1\leq l\leq
8,$ $\displaystyle s_{ijkl}\geq 0,\;(i,j,k)\in\Omega_{2},1\leq l\leq 8,$
$\displaystyle\lambda_{ijkl}\geq 0,\;(i,j,k)\in\Omega_{2},1\leq l\leq 8,$
$\displaystyle\sum_{1\leq i,j\leq 2n}G_{i,j}=0,$ $\displaystyle G\succeq 0.$

We use a tuning parameter $\rho$ for the trace in the objective function,
which can be used to balance obtaining a low-rank solution versus satisfying
the constraints. The tuning parameter $\rho$ can be chosen using cross-
validation or by selecting it so that the resulting solution has small
$(d+1)^{th}$ eigenvalue. As shown in section SM3 and section SM7, we observe
on synthetic and real data that the solution is robust to the choice of
$\rho$.

The theoretical results from Lemma 2.4 allow us to compute the distances
between homologous pairs from the pairwise distances $D_{ij}$. We recall that
we need to compute $\|v\|^{2}$ such that

$\det(T^{\prime}-8J\|v\|^{2})=0,$

where $T^{\prime}$ is an invertible matrix constructed from the pairwise
distance matrix by selecting a set of $2d+2$ indices. One step of computing
$\|v\|$ involves inverting $T^{\prime}$. Even if the error in the measurements
is small, noise can propagate and severely impact this computation. In order
to obtain a robust estimate of homolog-homolog distances, for each locus $i$,
we sample 100 $T^{\prime}$ matrices and obtain 100 solutions to the equation
for $\|v\|^{2}$. We then take the median of the solutions to be the homolog-
homolog distance for locus $i$ and use these homolog-homolog distances for the
evaluation of our algorithms on synthetic and real data in the following
section.

To solve the two convex optimization problems presented in this section for
the noiseless and noisy setting, we make use of the solver MOSEK implemented
in CVX within MATLAB. This results in the Gram matrix. In order to reconstruct
the coordinates of the genomic regions from the Gram matrix, we use an
eigenvector decomposition as also done in [46], namely: letting
$\gamma_{1},\ldots,\gamma_{d}$ be the top $d$ eigenvalues and
$\nu_{1},\ldots,\nu_{d}$ the corresponding eigenvectors of $G$, then

$x_{i}=(\sqrt{\gamma_{1}}\cdot\nu_{1,i},\ldots,\sqrt{\gamma_{d}}\cdot\nu_{d,i})\text{
and
}y_{i}=(\sqrt{\gamma_{1}}\cdot\nu_{1,n+i},\ldots,\sqrt{\gamma_{d}}\cdot\nu_{d,n+i})\text{
for }i=1,\ldots,n.$

Since we are interested in recovering the genome configuration in 3D, we use
$d=3$, thereby obtaining the desired 3D diploid configuration. We provide the
code for our algorithm at
https://github.com/uhlerlab/diploid-3D-reconstruction.

##  6 Evaluation on synthetic and real data

###  6.1 Synthetic data

We start by testing our method on simulated data. For this we construct three
different types of 3D structures: (a) a Brownian motion model using a standard
normal distribution to generate successive points; (b) points sampled
uniformly along a spiral with random translations sampled uniformly within
$(0,0.5)$ range and orientations sampled uniformly within
$(-\frac{\pi}{4},\frac{\pi}{4}$); (c) points sampled uniformly in a unit
sphere.

Performance of our method in the noiseless setting. For the 1D setting we
deduced in Section 2 that the pairwise distance constraints by themselves are
sufficient to identify the underlying 3D configuration. For the 2D setting we
proved in Section 3 that knowing additionally the distances between
neighboring beads leads to uniqueness. We here perform simulations in 3D since
this is the biologically relevant setting. These results are depicted in
Figure 5 with additional examples in Figure SM1. The input to our algorithm
are the pairwise distances (which are summed over homologs), all 3-way
distances, the distances between homologous loci, and the distances between
neighboring beads. In the noiseless setting considered here we solve the SDP
formulation in Equation 8. Figure 5 and Figure SM1 show that the true and
reconstructed structures highly overlap, thereby indicating that our
optimization formulation is able to recover the 3D structure of the full
diploid genome in the noiseless setting. When the 3-way distance constraints
are removed, the reconstructions are less aligned with the true structures.
This is shown in Figure 6, where we measure the root-mean-square deviation
(RMSD) between true and reconstructed 3D coordinates over $20$ trials. In line
with our theoretical results, these experimental results in the noiseless
setting indicate the importance of higher-order contact frequencies for
recovering the 3D diploid configuration, especially when the number of
chromosomes is larger.

(a)

(b)

(c)

Figure 5:  Examples of true and reconstructed points on simulated data. (a)
Brownian motion model. (b) Spirals. (c) Random points in a sphere. We generate
six chromosomes with in total of 120 domains, corresponding to three
homologous pairs with 20 domains per chromosome in the noiseless setting.
Solid lines / points correspond to true 3D coordinates and dashed lines /
unfilled points to reconstructions via our method. Each color represents a
different chromosome.

(a)

(b)

(c)

Figure 6: Performance of our method in the noiseless setting. Root-mean-square
deviation (RMSD) between true and reconstructed structure computed with and
without higher-order distance constraints. Simulated data was generated using
a Brownian motion model with (a) one (b) two and (c) three chromosomes. Mean
and standard deviation over $20$ trials are shown.

Performance of our method in the noisy setting. Next, we consider noisy
distance observations $D_{ij}=D_{ij}^{*}(1+\delta)$ and noisy 3-way distance
observations $D_{i_{1}i_{2}\ldots i_{k}}=D_{i_{1}i_{2}\ldots
i_{k}}^{*}(1+\delta)$ by sampling $\delta$ uniformly within
$(-\epsilon,\epsilon)$ as in [46], where $\epsilon$ is a given noise level.
For our simulations we sample a maximum of $1000$ 3-way distance constraints.
As shown in Figure SM2, we observe that the number of constraints does not
have a major effect on the reconstruction accuracy. While for all simulations
shown in this section, we set the tuning parameter $\rho=0.000001$, Figure SM3
shows that the performance is not significantly different when using different
choices of $\rho$.

In Figure 7 we numerically assess the accuracy of our predicted structure for
the Brownian motion model for different number of chromosomes (one, two, or
three) and different number of domains per chromosome ($10$ or $20$) by
computing the Spearman correlation between reconstructed and true pairwise
distances, similar to [46]. As expected, Figure 7 shows that when the noise
level increases, then the Spearman correlation between the original and
reconstructed configuration decreases. For the simulations with one
chromosome, the Spearman correlation is higher for $20$ domains than $10$.

(a)

(b)

(c)

Figure 7: Performance of our method in the noisy setting. Spearman correlation
under different noise levels for (a) one, (b) two and (c) three chromosomes.
Simulated data was generated using a Brownian motion model where each
chromosome has $10$ or $20$ domains. Mean and standard deviation over $20$
trials are shown.

###  6.2 Application to 3D diploid genome reconstruction

We apply our algorithm to the problem of reconstructing the diploid genome
from contact frequency data derived from experiments. We obtain pairwise and
3-way contact frequencies collected via SPRITE in human lymphoblastoid cells
from [33]. Since we aim to reconstruct the whole diploid genome, which
consists of approximately 6 billion base pairs, for computational reasons we
bin the contact frequencies in the SPRITE dataset into 10 Mega-base pair (Mb)
regions. While some previous studies considered higher resolutions, the
majority of the studies [4, 19, 36, 42, 46] did not attempt to reconstruct the
whole diploid genome and focused only on reconstructing one chromosome, thus
enabling them to consider higher resolutions.

After filtering out regions with a small number of total contacts, we obtain
514 unphased points on the chromosomes. We convert the pairwise contact
frequencies to pairwise distances using the previously observed relationship
$D_{ij}=F_{ij}^{-1/2}$ [36] and use Lemma 2.4 to obtain the distances between
homolog pairs from this data. As in our simulations in the noisy setting, we
randomly sample 1000 3-way distance constraints from all nonzero 3-way contact
frequencies (for the transformation from 3-way contact frequencies to 3-way
distances, see Section 4). Finally, we obtain the distances between
neighboring 10Mb beads by empirically evaluating the 3D reconstructions under
different input distances; see section SM4 and Figure SM4, Figure SM5.

Using the pairwise constraints, homolog-homolog constraints, neighboring bead
constraints, and 3-way distance constraints, we solve the SDP problem in
Equation 9 for the noisy setting and analyze the corresponding 3D coordinates.
Our diploid reconstruction is shown in Figure 8a. We compare this diploid
genome reconstruction to the 3D structure obtained via ChromSDE, shown in
Figure 8b obtained under the assumption that the observed contact frequencies
and the corresponding distances are a sum of four equal quantities, i.e.,
$\|x_{i}-x_{j}\|^{2},\|x_{i}-y_{j}\|^{2},\|y_{i}-x_{j}\|^{2}$, and
$\|y_{i}-y_{j}\|^{2}$ are equal. In Figure SM6, we show that the
reconstruction obtained using ChromSDE with equal distances does not
recapitulate known biology as described in the following paragraphs.

(a)

(b)

(c)

(d)

(e)

Figure 8: 3D Diploid Genome Reconstruction. Estimated 3D positions of all
chromosomes and their corresponding homologs at 10Mb resolution. 3D positions
obtained using (a) our method and (b) using ChromSDE with chromosomes colored
according to chromosome number. (c) Whole diploid organization obtained via
our method, colored by chromosome size. (d) Mean chromosome size as the
distance from the center increases. (e) The number of A compartments as the
distance from the center increases.

Experimental (imaging) studies have shown that chromosomes are organized by
size within the nucleus, with small chromosomes in the interior and larger
chromosomes on the periphery [3]. We colored each chromosome according to its
size and computed the mean chromosome size versus distance away from the
center. The results of the 3D configuration obtained using our method are
shown in Figures 8c and 8d and recapitulate prior studies: smaller chromosomes
are preferentially located in the center, whereas larger chromosomes are
preferentially on the periphery; see also section SM6, Figure SM7. This is
especially apparent for chromosomes 2 and 4, which are some of the largest
chromosomes, and in our reconstruction they are located on the periphery as
expected.

Experimental studies on the spatial organization of the genome have also shown
that the center of the nucleus is enriched in active compartments (known as A
compartments), while the periphery contains inactive compartments (known as B
compartments) [38]. From previously published data on the location of A and B
compartments along the genome in human lymphoblastoid cells [34], we counted
the number of A compartments per 10Mb bin. Then dividing our 3D reconstruction
into concentric circles of increasing radius away from the center, we found
the mean number of A compartments in each concentric circle. Figure 8e shows
that with increasing distance away from the center, the number of A
compartments decreases. Thus, our reconstruction recovers the experimentally
observed trend for A compartments to be preferentially located near the
nucleus center. As shown in Figures SM8 and SM9 in section SM7, we note that
our results are robust to the choice of the tuning parameter $\rho$ resulting
in biologically plausible configurations independent of the choice of $\rho$.

Currently, many studies such as [35] simply ignore the fact that the genome is
diploid and infer the 3D genome organization as if the data was collected from
a haploid organism, assuming that the homologous loci have the same 3D
structure. However, we show in section SM8, Figure SM10 that the haploid
distance matrices, computed by including only one copy of each of the
homologous loci, are different between the two copies with a mean Spearman
correlation of only 0.08. This shows that modeling the diploid aspect of the
genome provides valuable information regarding the 3D structure of each of the
homologs, which may be substantially different.

##  7 Discussion

In this paper, we proved that for a diploid organism the 3D genome structure
is not identifiable from pairwise distance measurements alone. This implies
that applying any algorithm for the reconstruction of the 3D genome structure
from typical chromosome conformation capture data for a diploid organism can
result in any of the infinitely many configurations with the same pairwise
contact frequencies. We showed that unique idenfiability is obtained using
distance constraints between neighboring genomic loci as well as 3-way
distance constraints in addition to the pairwise distance constraints that can
be obtained from typical contact frequency data. Distances between neighboring
genomic loci can be obtained empirically, e.g. from imaging studies, while
3-way distance constraints can be obtained from the most recently developed
sequencing-based methods for obtaining contact frequencies such as SPRITE
[33], C-walks [31] and GAM [2]. We also presented SDP formulations for
determining the 3D genome reconstruction both in the noiseless and the noisy
setting. Finally, we applied our algorithm to contact frequency data from
human lymphoblastoid cells collected using SPRITE and showed that our results
recapitulate known biological trends; in particular, in the 3D configuration
identified using our method, the small chromosomes are preferentially situated
in the interior of the cell nucleus, while the larger chromosomes are
preferentially situated at the periphery of the cell nucleus. In addition, in
the 3D configuration identified using our method the number of A domains is
higher in the interior versus the periphery, which is in line with
experimental results. Our work shows the importance of higher-order contact
frequencies that can be measured using SPRITE [33], C-walks [31] and GAM [2]
for obtaining the 3D organization of the genome in diploid organisms. This is
particularly relevant for the reconstruction of cancer genomes, where copy
number variations are frequent and hence the genome may contain even more than
two copies of each locus.

We conjecture that identifiability of the 3D genome structure can also be
achieved by replacing the higher-order contact constraints by distance
constraints to the center of the cell nucleus. Such constraints are also
biologically relevant, since these distances can be measured via imaging
experiments, or inferred by measuring whether a particular locus is in a
lamin-associated domain or a telomere, both of which tend to lie at the
boundary of the cell nucleus [10, 18, 41]. Another future research direction
is the development of specialized solvers to enable reconstruction of the
genome at higher resolution. In this study we used a 10Mbp resolution due to
the computational constraints imposed by SDP solvers. Finally, the theoretical
results in this paper build on the assumption that distances are inverses of
square roots of pairwise and higher-order contact frequencies. An interesting
future research direction is to develop a method for estimating the map
between higher-order contact frequencies and distances, and then prove
identifiability as well as build reconstruction algorithms for these different
maps.

## Acknowldegements

We thank Mohab Safey El Din for helpful discussions.

## References

  * [1] A. Y. Alfakih, A. Khandani, and H. Wolkowicz, Solving Euclidean distance matrix completion problems via semidefinite programming, Computational Optimization and Applications, 12 (1999), pp. 13–30.
  * [2] R. A. Beagrie, A. Scialdone, M. Schueler, D. C. A. Kraemer, M. Chotalia, S. Q. Xie, M. Barbieri, I. de Santiago, L.-M. Lavitas, M. R. Branco, et al., Complex multi-enhancer contacts captured by genome architecture mapping, Nature, 543 (2017), p. 519.
  * [3] A. Bolzer, G. Kreth, I. Solovei, D. Koehler, K. Saracoglu, C. Fauth, S. Müller, R. Eils, C. Cremer, M. R. Speicher, et al., Three-dimensional maps of all chromosomes in human male fibroblast nuclei and prometaphase rosettes, PLoS Biology, 3 (2005), p. e157.
  * [4] A. G. Cauer, G. Yardimci, J.-P. Vert, N. Varoquaux, and W. S. Noble, Inferring diploid 3D chromatin structures from Hi-C data, bioRxiv:644294, (2019).
  * [5] L. Cayton and S. Dasgupta, Robust Euclidean embedding, in Proceedings of the 23rd International Conference on Machine Learning, ICML ’06, New York, NY, USA, 2006, ACM, pp. 169–176.
  * [6] V. Chandrasekaran, B. Recht, P. A. Parrilo, and A. S. Willsky, The convex geometry of linear inverse problems, Foundations of Computational Mathematics, 12 (2012), pp. 805–849.
  * [7] . G. P. Consortium et al., An integrated map of genetic variation from 1,092 human genomes, Nature, 491 (2012), pp. 56–65.
  * [8] . G. P. Consortium et al., A global reference for human genetic variation, Nature, 526 (2015), pp. 68–74.
  * [9] T. F. Cox and M. A. A. Cox, Multidimensional scaling, Chapman and Hall / CRC, 2000.
  * [10] L. Crabbe, A. J. Cesare, J. M. Kasuboski, J. A. J. Fitzpatrick, and J. Karlseder, Human telomeres are tethered to the nuclear envelope during postmitotic nuclear assembly, Cell Reports, 2 (2012), pp. 1521–1529.
  * [11] J. Dekker, Gene regulation in the third dimension, Science, 319 (2008), pp. 1793–1794.
  * [12] J. Dekker, K. Rippe, M. Dekker, and N. Kleckner, Capturing chromosome conformation, Science, 295 (2002), pp. 1306–1311.
  * [13] M. Di Pierro, B. Zhang, E. L. Aiden, P. G. Wolynes, and J. N. Onuchic, Transferable model for chromosome architecture, Proceedings of the National Academy of Sciences, 113 (2016), pp. 12168–12173.
  * [14] Z. Duan, M. Andronescu, K. Schutz, S. McIlwain, Y. J. Kim, C. Lee, J. Shendure, S. Fields, C. A. Blau, and W. S. Noble, A three-dimensional model of the yeast genome, Nature, 465 (2010), pp. 363–367, http://dx.doi.org/10.1038/nature08973.
  * [15] H. Fang and D. P. O’Leary, Euclidean distance matrix completion problems, Optimization Methods and Software, 27 (2012), pp. 695–717.
  * [16] M. Fazel, H. Hindi, and S. P. Boyd, Log-det heuristic for matrix rank minimization with applications to Hankel and Euclidean distance matrices, in Proceedings of the 2003 American Control Conference, vol. 3, IEEE, 2003, pp. 2156–2162.
  * [17] M. Fazel, H. Hindi, S. P. Boyd, et al., A rank minimization heuristic with application to minimum order system approximation, in Proceedings of the American Control Conference, vol. 6, Citeseer, 2001, pp. 4734–4739.
  * [18] L. Guelen, L. Pagie, E. Brasset, W. Meuleman, M. B. Faza, W. Talhout, B. H. Eussen, A. de Klein, L. Wessels, W. de Laat, et al., Domain organization of human chromosomes revealed by mapping of nuclear lamina interactions, Nature, 453 (2008), pp. 948–951.
  * [19] M. Hu, K. Deng, Z. Qin, J. Dixon, S. Selvaraj, J. Fang, B. Ren, and J. S. Liu, Bayesian inference of spatial organizations of chromosomes, PLoS Computational Biology, 9 (2013), p. e1002893.
  * [20] J. R. Hughes, N. Roberts, S. McGowan, D. Hay, E. Giannoulatou, M. Lynch, M. De Gobbi, S. Taylor, R. Gibbons, and D. R. Higgs, Analysis of hundreds of cis-regulatory landscapes at high resolution in a single, high-throughput experiment, Nature Genetics, 46 (2014), pp. 205–212.
  * [21] R. Jungmann, M. S. Avendaño, J. B. Woehrstein, M. Dai, W. M. Shih, and P. Yin, Multiplexed 3d cellular super-resolution imaging with dna-paint and exchange-paint, Nature Methods, 11 (2014), p. 313.
  * [22] N. Krislock, Semidefinite facial reduction for low-rank Euclidean distance matrix completion, PhD thesis, University of Waterloo, 2010, http://hdl.handle.net/10012/5093.
  * [23] A. Lesne, J. Riposo, P. Roger, A. Cournac, and J. Mozziconacci, 3D genome reconstruction from chromosomal contacts, Nature Methods, 11 (2014), pp. 1141–1143.
  * [24] E. Lieberman-Aiden, N. L. v. Berkum, L. Williams, M. Imakaev, T. Ragoczy, A. Telling, I. Amit, B. R. Lajoie, P. J. Sabo, M. O. Dorschner, R. Sandstrom, B. Bernstein, M. A. Bender, M. Groudine, A. Gnirke, J. Stamatoyannopoulos, L. A. Mirny, E. S. Lander, and J. Dekker, Comprehensive mapping of long-range interactions reveals folding principles of the human genome, Science, 326 (2009).
  * [25] E. Lieberman-Aiden, N. L. Van Berkum, L. Williams, M. Imakaev, T. Ragoczy, A. Telling, I. Amit, B. R. Lajoie, P. J. Sabo, M. O. Dorschner, R. Sandstrom, B. Bernstein, M. A. Bender, M. Groudine, A. Gnirke, J. Stamatoyannopoulos, L. A. Mirny, E. S. Lander, and J. Dekker, Comprehensive mapping of long-range interactions reveals folding principles of the human genome, science, 326 (2009), pp. 289–293.
  * [26] F. Lu, S. Keleş, S. J. Wright, and G. Wahba, Framework for kernel regularization with application to protein clustering, Proceedings of the National Academy of Sciences, 102 (2005), pp. 12332–12337.
  * [27] L. A. Mirny, The fractal globule as a model of chromatin architecture in the cell, Chromosome research, 19 (2011), pp. 37–51.
  * [28] B. Mishra, G. Meyer, and R. Sepulchre, Low-rank optimization for distance matrix completion, in 2011 50th IEEE Conference on Decision and Control and European Control Conference, 2011, pp. 4455–4460.
  * [29] I. Müller, S. Boyle, R. H. Singer, W. A. Bickmore, and J. R. Chubb, Stable morphology, but dynamic internal reorganisation, of interphase human chromosomes in living cells, PloS One, 5 (2010), p. e11560.
  * [30] G. Nir, I. Farabella, C. P. Estrada, C. G. Ebeling, B. J. Beliveau, H. M. Sasaki, S. H. Lee, S. C. Nguyen, R. B. McCole, S. Chattoraj, et al., Walking along chromosomes with super-resolution imaging, contact maps, and integrative modeling, PLoS Genetics, 14 (2018), p. e1007872.
  * [31] P. Olivares-Chauvet, Z. Mukamel, A. Lifshitz, O. Schwartzman, N. O. Elkayam, Y. Lubling, G. Deikus, R. P. Sebra, and A. Tanay, Capturing pairwise and multi-way chromosomal conformations using chromosomal walks, Nature, 540 (2016), pp. 296–300.
  * [32] Y. Qi and B. Zhang, Predicting three-dimensional genome organization with chromatin states, PLoS computational biology, 15 (2019), p. e1007024.
  * [33] S. A. Quinodoz, N. Ollikainen, B. Tabak, A. Palla, J. M. Schmidt, E. Detmar, M. M. Lai, A. A. Shishkin, P. Bhat, Y. Takei, et al., Higher-order inter-chromosomal hubs shape 3D genome organization in the nucleus, Cell, 174 (2018), pp. 744–757.
  * [34] S. S. P. Rao, M. H. Huntley, N. C. Durand, E. K. Stamenova, I. D. Bochkov, J. T. Robinson, A. L. Sanborn, I. Machol, A. D. Omer, E. S. Lander, et al., A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping, Cell, 159 (2014), pp. 1665–1680.
  * [35] L. Rieber and S. Mahony, miniMDS: 3D structural inference from high-resolution Hi-C data, Bioinformatics, 33 (2017), pp. i261–i266.
  * [36] M. Rousseau, J. Fraser, M. A. Ferraiuolo, J. Dostie, and M. Blanchette, Three-dimensional modeling of chromatin structure from interaction frequency data using Markov chain Monte Carlo sampling, BMC Bioinformatics, 12 (2011), p. 414.
  * [37] M. Simonis, P. Klous, E. Splinter, Y. Moshkin, R. Willemsen, E. de Wit, B. van Steensel, and W. de Laat, Nuclear organization of active and inactive chromatin domains uncovered by chromosome conformation capture–on-chip (4C), Nature Genetics, 38 (2006), pp. 1348–1354.
  * [38] T. J. Stevens, D. Lando, S. Basu, L. P. Atkinson, Y. Cao, S. F. Lee, M. Leeb, K. J. Wohlfahrt, W. Boucher, A. O’Shaughnessy-Kirwan, et al., 3D structures of individual mammalian genomes studied by single-cell Hi-C, Nature, 544 (2017), p. 59.
  * [39] C. Uhler and G. V. Shivashankar, Chromosome intermingling: mechanical hotspots for genome regulation, Trends in Cell Biology, 27 (2017), pp. 810–819.
  * [40] C. Uhler and G. V. Shivashankar, The regulation of genome organization and gene expression by nuclear mechanotransduction, Nature Reviews Molecular Cell Biology, 18 (2017), pp. 717–727.
  * [41] B. Van Steensel and A. S. Belmont, Lamina-associated domains: links with chromosome architecture, heterochromatin, and gene repression, Cell, 169 (2017), pp. 780–791.
  * [42] N. Varoquaux, F. Ay, W. S. Noble, and J.-P. Vert, A statistical approach for inferring the 3D structure of the genome, Bioinformatics, 30 (2014), pp. i26–i33.
  * [43] H. Wang, X. Xu, C. M. Nguyen, Y. Liu, Y. Gao, X. Lin, T. Daley, N. H. Kipniss, M. La Russa, and L. S. Qi, CRISPR-mediated programmable 3D genome positioning and nuclear organization, Cell, 175 (2018), pp. 1405–1417.
  * [44] K. Q. Weinberger, F. Sha, Q. Zhu, and L. K. Saul, Graph Laplacian regularization for large-scale semidefinite programming, in Advances in Neural Information Processing Systems, 2007, pp. 1489–1496.
  * [45] L. Zhang, G. Wahba, and M. Yuan, Distance shrinkage and Euclidean embedding via regularized kernel estimation, Journal of the Royal Statistical Society: Series B, 78 (2016), pp. 849–867.
  * [46] Z. Zhang, G. Li, K.-C. Toh, and W.-K. Sung, 3D Chromosome Modeling with Semi-Definite Programming and Hi-C Data, Journal of Computational Biology, 20 (2013).