Fundamental Limits of Multiple Sequence Reconstruction from Substrings

The problem of reconstructing a sequence from the set of its length-$k$ substrings (or $k$-mers) has a long history, dating back to the development of the sequencing-by-hybridization (SBH) technology in the late 1980s [1]. Due to this technological advance, significant attention was originally devoted to understanding when the set of $k$-mers uniquely determines the underlying sequence, and to designing algorithms to perform the reconstruction efficiently [2, 3, 4, 5].

Inspired by Ukkonen’s 1992 paper [2], which characterized how large $k$ needs to be to allow reconstruction, a more recent line of work considered this problem from an information-theoretic perspective [6, 7, 8, 9, 10]. In particular, [7] considered a random source ${\bf x}$ (modeling a genome sequence), and provided conditions on how large $k$ needs to be (and how many length-$k$ substrings need to be observed) in order for the reconstruction of a source sequence ${\bf x}$ to be information-theoretically feasible.

In this work, we consider a multiple-source version of this problem. We assume that there are $m$ random source sequences ${\bf x}_{1},\dots,{\bf x}_{m}$, each of length $n$, and we wish to reconstruct them from the union of their sets of $k$-mers. This multiple-source setting can be motivated by modern sequencing assays such as metagenomics and RNA-seq, where one observes substrings from a mixture of distinct source sequences (different microbial genomes in the case of metagenomics and different mRNA molecules in the case of RNA-seq). It can also be motivated by the fact that many bioinformatics algorithms first convert a set of reads into the set of their constituent $k$-mers and then operate in $k$-mer space. It is thus relevant to understand when such sets fully preserve the information in the original reads. Our setting is also connected with the problem of DNA-based data storage, where one wishes to recover a set of $m$ synthesized DNA molecules by jointly sequencing them. Our work can thus be seen as an uncoded version of the recent work on multi-strand reconstruction from substrings [11], or a random coding approach to it (see also “Related work” below).

Our goal is to characterize how large $k$ needs to be to allow perfect reconstruction of all $m$ source sequences. We assume each source is i.i.d. ${\rm Bern}(\tfrac{1}{2})$ and consider an asymptotic regime where $m=n^{\alpha}$ and $k=\beta\log n$, for some constants $\alpha,\beta>0$. We aim to characterize for which $(\alpha,\beta)$ pairs the set of $k$-mers uniquely determines the set of source sequences $\{{\bf x}_{1},\dots,{\bf x}_{m}\}$ with vanishing error probability. Our main result, illustrated in Figure 1, characterizes the $(\alpha,\beta)$ feasibility region almost completely. We establish this feasibility region by carefully analyzing the (random) de Bruijn graph formed by the $k$-mer set. In this graph, the set of source sequences corresponds to a set of paths that cover all edges. Our analysis shows that there are feasible $(\alpha,\beta)$ pairs where the $k$-mer de Bruijn graph is fairly complex, due to the presence of repeats across the source sequences, but nevertheless there is a unique way to decompose the graph into $m$ source paths.

A multi-sequence version of the coding problem in [12] was considered in [13]. Similar to our paper, they provide conditions on how $k$ needs to scale with $n$ and $m$ to guarantee that the asymptotic rate of all $k$-mer-reconstructible codes approaches 1. Our work can be seen as an uncoded counterpart to [13], where we characterize how $k$ needs to scale with $n$ and $m$ to allow the reconstruction of random source sequences with vanishing failure probability. As an important distinction, we notice that the code constructions in [13] guarantee that there are no repeats across the $m$ sequences to be reconstructed, while in our setting, $m$ source sequences with repeats may still be reconstructible for some $(\alpha,\beta)$.

We consider $m$ source sequences ${\bf x}_{1},\dots,{\bf x}_{m}$, each of which is generated as an i.i.d. ${\rm Ber}(\tfrac{1}{2})$ length-$n$ sequence. The goal is to reconstruct the set of source sequences ${\mathcal{X}}:=\{{\bf x}_{1},\dots,{\bf x}_{m}\}$ from the set of all of their $(k+1)$-mers, i.e.,

$\displaystyle{\mathcal{Y}}({\mathcal{X}}):=\{{\bf x}_{i}[j:j+k]:1\leq i\leq m,\,1\leq j\leq n-k\},$

where ${\bf x}_{i}[a:b]$ is the substring of ${\bf x}_{i}$ between (and including) the $a$th and $b$th symbols. Here, we use $(k+1)$-mers rather than $k$-mers for notational convenience, and the final results are not affected. When clear from context, we write ${\mathcal{Y}}$ for ${\mathcal{Y}}({\mathcal{X}})$. We will also use the notation ${\bf x}_{i}(a)$ to refer to the $a$th $k$-mer of ${\bf x}_{i}$, for $1\leq a\leq n-k+1$.

We consider an asymptotic regime where $m=n^{\alpha}$ and $k=\beta\log n$, for positive constants $\alpha$ and $\beta$. As it turns out, this is the asymptotic regime of interest, where the fundamental limits are non-trivial. We say that ${\mathcal{Y}}$ reconstructs ${\mathcal{X}}$ if ${\mathcal{X}}$ is the unique set of $m$ length-$n$ sequences (up to relabeling) that could have generated ${\mathcal{Y}}$. It is then natural to define a notion of feasibility as follows.

The goal is to characterize the feasibility region ${\mathcal{F}}\subset{\mathbb{R}}_{+}^{2}$, defined as the set of all feasible $(\alpha,\beta)$ pairs.

An equivalent and useful representation of the $k$-mer set ${\mathcal{Y}}$ is its de Bruijn graph $G({\mathcal{Y}})$. The de Bruijn graph $G({\mathcal{Y}})$ can be obtained by taking the set of $k$-mers in ${\mathcal{Y}}$ (i.e., length-$k$ prefixes and suffixes of the $(k+1)$-mers in ${\mathcal{Y}}$) as the node set, and connecting two $k$-mers with a directed edge if they are the prefix and suffix of a $(k+1)$-mer in ${\mathcal{Y}}$. Notice that the true set of sequences ${\mathcal{X}}$ can be seen as a set of paths on $G({\mathcal{Y}})$ that cover all the edges (edges can be used multiple times).

Our main result is to characterize, almost completely, the feasibility region ${\mathcal{F}}$ as follows:

As illustrated in Figure 1, Theorem 1 nearly completely characterizes ${\mathcal{F}}$. In particular, for $\alpha>1$, $(\alpha,\beta)\in{\mathcal{F}}$ if $\beta>2\alpha+1$, and $(\alpha,\beta)\notin{\mathcal{F}}$ if $\beta<2\alpha+1$. For $\alpha<1$ there is a small uncharacterized region of $(\alpha,\beta)$ pairs.

The rest of the paper is organized as follows. In Section III, as a warm-up, we characterize the repeat-free region, a (strict) subset of ${\mathcal{F}}$ where there are no length-$k$ repeats across the source sequences (and the source strings are thus reconstructible). In Section IV, we prove the inner bound part of Theorem 1, and in Section V, we prove the outer bound part of Theorem 1.

We start by considering a simple scenario where $\mathcal{X}$ can be recovered correctly: when each $k$-mer is unique across all source sequences. In this case, it is well known that the de Bruijn graph $G({\mathcal{Y}})$ will contain $m$ disjoint paths, each of which corresponds to one of the ${\bf x}_{i}$s, and ${\mathcal{X}}$ is thus reconstructible. We prove the following result.

This region is marked in Figure 1. Notice that it does not include all pairs of feasible $(\alpha,\beta)$ pairs.

To prove Lemma 1, we note that, if $i\neq j$,

$\displaystyle\Pr\left({\bf x}_{i}(a)={\bf x}_{j}(b)\right)=2^{-k}=2^{-\beta\log n}=n^{-\beta}.$

Hence, the probability that any two strings ${\bf x}_{i}$ and ${\bf x}_{j}$, for $i\neq j$ have at least one $k$-mer in common is upper-bounded using the union bound by

$\displaystyle\binom{m}{2}n^{2}n^{-\beta}\leq m^{2}n^{2-\beta}\leq n^{2\alpha+2-\beta}.$

Additionally, we need to consider the case where a single string contains a repeat $k$-mer. In this case, not all of the $k$-mers are independent due to potential overlaps. However, a careful analysis reveals that, even if $|b-a|<k$ and ${\bf x}_{i}(a)$ and ${\bf x}_{i}(b)$ overlap, we still have $\Pr\left({\bf x}_{i}(a)={\bf x}_{i}(b)\right)=n^{-\beta}$ (see Section A for a similar analysis). We conclude that the probability that there is repeated $k$-mer within some ${\bf x}_{i}$ is at most $mn^{2}n^{-\beta}=n^{\alpha+2-\beta}$, and the probability that any $k$-mer is repeated across or within ${\bf x}_{1},\dots,{\bf x}_{m}$ is upper-bounded by

$\displaystyle n^{2\alpha+2-\beta}+n^{\alpha+2-\beta},$

which tends to zero for $\beta>2\alpha+2$. A converse result can also be shown; i.e., for $\beta<2\alpha+2$, there will be repeats with probability tending to $1$ (but we do not describe that). In the next section, we show that there are feasible $(\alpha,\beta)$ pairs that are in this repeat-abundant regime.

In order to characterize ${\mathcal{F}}$, we define the error event ${\mathcal{E}}=\{{\mathcal{Y}}\text{ does not reconstruct }{\mathcal{X}}\}$, and seek to characterize the asymptotic behavior of $\Pr({\mathcal{E}})$ for a given choice of $(\alpha,\beta)$. To do so, we will need to define several events related to repeats:

	$\displaystyle{\mathcal{A}}=\{$	$\displaystyle\exists i,a,b:{\bf x}_{i}(a)={\bf x}_{i}(b)\}$
	$\displaystyle{\mathcal{B}}=\{$	$\displaystyle\exists i,j,l,a,b,c,d:{\bf x}_{i}(a)={\bf x}_{j}(c),{\bf x}_{i}(b)={\bf x}_{l}(d),$
		$\displaystyle 0\leq b-a<k,\text{ and either }\{j\neq l\}\text{ or }$
		$\displaystyle\{j=l\text{ and }d-c\neq b-a\}\}$
	$\displaystyle{\mathcal{C}}=\{$	$\displaystyle\exists i,j,a:{\bf x}_{i}(1)={\bf x}_{j}(a)\}\cup\{\exists i,j,a:{\bf x}_{i}(n^{\prime})={\bf x}_{j}(a)\}$
	$\displaystyle{\mathcal{D}}=\{$	$\displaystyle\exists i,j,a:{\bf x}_{i}(a)={\bf x}_{j}(a)\}$

where we let $n^{\prime}=n-k+1$ for conciseness.

Event ${\mathcal{A}}$ is the event that there is at least one intra-sequence repeat. Event ${\mathcal{B}}$ is the event that there are two pairs of repeats with an overlap (i.e., ${\bf x}_{i}(a)$ and ${\bf x}_{i}(b)$ overlap in ${\bf x}_{i}$), but they are not simply a longer length-$t$ repeat for $t>k$. Notice that ${\mathcal{B}}$ includes the event that a triple repeat ${\bf x}_{i}(a)={\bf x}_{j}(b)={\bf x}_{l}(c)$ occurs. Event ${\mathcal{C}}$ is the event that the first or last $k$-mer of some ${\bf x}_{i}$ appears somewhere else. Event ${\mathcal{D}}$ is the event that the $a$th $k$-mer of two sequences ${\bf x}_{i}$ and ${\bf x}_{j}$ is the same. Notice that events ${\mathcal{A}}$, ${\mathcal{B}}$ and ${\mathcal{C}}$ are not error events per se (they each do not imply ${\mathcal{E}}$), but they will be useful in the analysis. On the other hand, ${\mathcal{D}}\Rightarrow{\mathcal{E}}$, since given ${\mathcal{D}}$ it is possible to replace ${\bf x}_{i}$ and ${\bf x}_{j}$ with two incorrect sequences

	$\displaystyle\tilde{\bf x}_{i}={\bf x}_{i}[1:a-1]{\bf x}_{j}[a:n]$
	$\displaystyle\tilde{\bf x}_{j}={\bf x}_{j}[1:a-1]{\bf x}_{i}[a:n].$

We will upper bound the error probability as

	$\displaystyle\Pr({\mathcal{E}})$	$\displaystyle\leq\Pr({\mathcal{A}}\cup{\mathcal{B}}\cup{\mathcal{C}}\cup{\mathcal{D}}\cup{\mathcal{E}})$
		$\displaystyle\leq\Pr({\mathcal{A}})+\Pr({\mathcal{B}})+\Pr({\mathcal{C}})+\Pr({\mathcal{D}})+$
		$\displaystyle\quad+\Pr({\mathcal{E}}\cap\bar{\mathcal{A}}\cap\bar{\mathcal{B}}\cap\bar{\mathcal{C}}\cap\bar{\mathcal{D}})$		(1)

and analyze the conditions for each term to go to zero. First we notice that, by the union bound,

$\displaystyle\Pr({\mathcal{A}})\leq mn^{2}n^{-\beta}=n^{\alpha+2-\beta},$

(2)

from which we see that $\Pr({\mathcal{A}})\to 0$ if $\beta>\alpha+2$. Similarly,

	$\displaystyle\Pr({\mathcal{B}})\leq m^{3}n^{3}2kn^{-2\beta}+m^{2}n^{2}(2k)^{2}n^{-2\beta}$
	$\displaystyle\quad\quad\quad\leq n^{3\alpha+3-2\beta+o(1)},$		(3)
	$\displaystyle\Pr({\mathcal{C}})\leq 2m^{2}nn^{-\beta}=2n^{2\alpha+1-\beta},$		(4)
	$\displaystyle\Pr({\mathcal{D}})\leq m^{2}nn^{-\beta}=4n^{2\alpha+1-\beta}.$		(5)

We point out that to establish (3), we need a careful analysis of the case where ${\bf x}_{j}(c)$ and ${\bf x}_{j}(d)$ also overlap, presented in Lemma 4.

As a result we see that, if $\beta>\max(\alpha+2,\tfrac{3}{2}\alpha+\tfrac{3}{2},2\alpha+1)$, $\Pr({\mathcal{A}})+\Pr({\mathcal{B}})+\Pr({\mathcal{C}})+\Pr({\mathcal{D}})\to 0$. What remains is to bound $\Pr({\mathcal{E}}\cap\bar{\mathcal{A}}\cap\bar{\mathcal{B}}\cap\bar{\mathcal{C}}\cap\bar{\mathcal{D}})$.

Now suppose $\bar{\mathcal{A}}\cap\bar{\mathcal{B}}\cap\bar{\mathcal{C}}\cap\bar{\mathcal{D}}$ holds. Let the multiplicity $\mu_{{\mathcal{X}}}({\bf v})$ of an node ${\bf v}$ in $G({\mathcal{Y}})$ be the number of times it is traversed by paths in ${\mathcal{X}}$. Notice that $\bar{\mathcal{B}}$ guarantees that no node can be traversed three or more times, so given $\bar{\mathcal{B}}$, node multiplicities can only be $1$ or $2$.

The proof of the following lemma is in Appendix B.

Now suppose ${\mathcal{E}}$ occurs. By definition, there is at least one $\tilde{\mathcal{X}}\neq{\mathcal{X}}$ such that ${\mathcal{Y}}(\tilde{\mathcal{X}})={\mathcal{Y}}({\mathcal{X}})$. Let $\tilde{\mathcal{X}}$ be one such choice with minimum set difference $|{\mathcal{X}}-\tilde{\mathcal{X}}|$. Let $|{\mathcal{X}}-\tilde{\mathcal{X}}|=c$. By Lemma 2, we know that for all nodes ${\bf v}$,

$\displaystyle\mu_{{\mathcal{X}}}({\bf v})=\mu_{\tilde{\mathcal{X}}}({\bf v}):=\mu_{G}({\bf v}).$

(6)

Additionally,

	$\displaystyle\mu_{{\mathcal{X}}\cap\tilde{{\mathcal{X}}}}({\bf v})+\mu_{{\mathcal{X}}-\tilde{{\mathcal{X}}}}({\bf v})=\mu_{{\mathcal{X}}}({\bf v})$	$\displaystyle=\mu_{G}({\bf v})$		(7)
	$\displaystyle\mu_{{\mathcal{X}}\cap\tilde{{\mathcal{X}}}}({\bf v})+\mu_{\tilde{{\mathcal{X}}}-{\mathcal{X}}}({\bf v})=\mu_{\tilde{\mathcal{X}}}({\bf v})$	$\displaystyle=\mu_{G}({\bf v}),$		(8)

implying that $\mu_{{\mathcal{X}}-\tilde{{\mathcal{X}}}}({\bf v})=\mu_{\tilde{{\mathcal{X}}}-{\mathcal{X}}}({\bf v})$.

We now build a “difference de Bruijn graph” $\tilde{G}({\mathcal{Y}},{\mathcal{X}},\tilde{\mathcal{X}})$ on the $(k+1)$-mers in ${\mathcal{Y}}$ that are present in ${\mathcal{X}}-\tilde{\mathcal{X}}$. It is easy to see that $\tilde{G}({\mathcal{Y}},{\mathcal{X}},\tilde{\mathcal{X}})$ is the same as the subgraph of $G({\mathcal{Y}})$ induced by the edges corresponding to $(k+1)$-mers in ${\mathcal{X}}-\tilde{\mathcal{X}}$ with nodes $\{{\bf v}:\mu_{{\mathcal{X}}-\tilde{\mathcal{X}}}({\bf v})>0\}$. Note that the nodes and their multiplicities are equivalent to those of the set $\{{\bf v}:\mu_{\tilde{\mathcal{X}}-{\mathcal{X}}}({\bf v})>0\}$, and the set of edges in $\tilde{G}({\mathcal{Y}},{\mathcal{X}},\tilde{\mathcal{X}})$ must also correspond to all of the $(k+1)$-mers in ${\mathcal{Y}}$ that are present in $\tilde{\mathcal{X}}-{\mathcal{X}}$. Therefore $\tilde{G}({\mathcal{Y}},{\mathcal{X}},\tilde{\mathcal{X}})$ is the same as the de Bruijn graph on the $(k+1)$-mers in ${\mathcal{Y}}$ that are present in $\tilde{\mathcal{X}}-{\mathcal{X}}$; i.e., $\tilde{G}({\mathcal{Y}},{\mathcal{X}},\tilde{\mathcal{X}})=\tilde{G}({\mathcal{Y}},\tilde{\mathcal{X}},{\mathcal{X}})$. Notice that $\tilde{G}({\mathcal{Y}},{\mathcal{X}},\tilde{\mathcal{X}})$ is a de Bruijn graph for both the set of $c$ missing true paths ${\mathcal{X}}-\tilde{\mathcal{X}}=\{{\bf x}_{1},\dots,{\bf x}_{c}\}$ and for the set of incorrectly reconstructed paths $\tilde{\mathcal{X}}-{\mathcal{X}}=\{\tilde{\bf x}_{1},\dots,\tilde{\bf x}_{c}\}$.

Define a maximal shared subpath in $\tilde{G}({\mathcal{Y}},{\mathcal{X}},\tilde{\mathcal{X}})$ as a directed sequence $\mathbf{p}$ of nodes with maximal length such that $\mu_{{\mathcal{X}}-\tilde{\mathcal{X}}}({\bf v})=2$ for all nodes ${\bf v}\in\mathbf{p}$. We have the following lemma.

Figure 2 shows examples of the graph $\tilde{G}({\mathcal{Y}},{\mathcal{X}},\tilde{\mathcal{X}})$ given $\bar{\mathcal{A}},\bar{\mathcal{B}},\bar{\mathcal{C}},\bar{\mathcal{D}}$, for $c=2$ and $c=4$. Notice that we have $c$ paths and $c$ maximal shared subpaths. Also notice that, in both examples, there exist two choices for the set of $c$ paths.

Lemma 3 allows us to analyze $\Pr({\mathcal{E}})$ by partitioning it into events ${\mathcal{E}}_{c}:={\mathcal{E}}\cap\{|\tilde{\mathcal{X}}-{\mathcal{X}}|=c\}$. Notice that a maximal shared subpath between ${\bf x}_{i}$ and ${\bf x}_{j}$ occurs with probability at most $n^{-\beta}$. Moreover, by the minimality of $\tilde{\mathcal{X}}$, each $\tilde{\bf x}_{i}$ must contain at least one maximal shared subpath. By the union bound, we conclude that

	$\displaystyle\Pr($	$\displaystyle{\mathcal{E}}_{c}\cap\bar{\mathcal{A}}\cap\bar{\mathcal{B}}\cap\bar{\mathcal{C}}\cap\bar{\mathcal{D}})$
		$\displaystyle\leq\binom{m}{c}n^{c}(cn)^{c}n^{-\beta c}\leq m^{c}n^{2c}n^{-\beta c}=n^{c(\alpha+2-\beta)},$

where we used the fact that $\tbinom{m}{c}\leq(m/c)^{c}$. Summing over all values of $c$, we obtain

	$\displaystyle\Pr({\mathcal{E}}\cap\bar{\mathcal{A}}\cap\bar{\mathcal{B}}\cap\bar{\mathcal{C}}\cap\bar{\mathcal{D}})$	$\displaystyle\leq\sum_{c=2}^{m}n^{c(\alpha+2-\beta)}$
		$\displaystyle=\frac{n^{2(\alpha+2-\beta)}(1-n^{(m-1)(\alpha+2-\beta)})}{1-n^{\alpha+2-\beta}},$

which tends to zero as $n\to\infty$ provided that $\beta>\alpha+2$. Plugging this back into (1), we conclude that any $(\alpha,\beta)$ pair with $\beta>\max(\alpha+2,\tfrac{3}{2}\alpha+\tfrac{3}{2},2\alpha+1)=\max(\alpha+2,2\alpha+1)$ is feasible. In other words, we showed that

$\displaystyle\{(\alpha,\beta):\beta>\max(\alpha+2,2\alpha+1)\}\subseteq{\mathcal{F}},$

concluding the proof of the achievability part of Theorem 1.

In this section, we will prove the outer bound part of Theorem 1. First, we will show that $\Pr({\mathcal{E}})\not\to 0$ as $n\to\infty$, if $\beta<2\alpha+1$. Second, we will show that $\Pr({\mathcal{E}})\not\to 0$ if $\beta<\alpha+\tfrac{3}{2}$. To do so, we will show that in those regimes, with positive probability, there will be multiple reconstructions.

Define $Z_{i,j,a}=\mathbbm{1}_{{\bf x}_{i}(a)={\bf x}_{j}(a)}$ for $i<j$; i.e., $Z_{i,j,a}$ indicates whether ${\bf x}_{i}$ and ${\bf x}_{j}$ share a $k$-mer at a common position $a$. Let $V=\sum_{i=1}^{m-1}\sum_{j=i+1}^{m}\sum_{a=1}^{n-k+1}Z_{i,j,a}$, and note that ${\mathcal{D}}=\{V>0\}$. In particular, since any $Z_{i,j,a}=1$ implies that there are multiple valid reconstructions, we have that $\Pr({\mathcal{E}})\geq\Pr(V>0)$. By the Paley-Zygmund inequality, we then have

$\displaystyle\Pr({\mathcal{E}})\geq\Pr(V>0)\geq\frac{\mathbb{E}[V]^{2}}{\mathbb{E}[V^{2}]}.$

(9)

It is straightforward to see $\mathbb{E}[Z_{i,j,a}]=n^{-\beta}$, and by the linearity of expectation,

$$\mathbb{E}[V]=\binom{m}{2}(n-k+1)n^{-\beta}.$$

The computation of $E[V^{2}]$ requires more care in handling $\Pr({\bf x}_{i}(a)={\bf x}_{j}(a),{\bf x}_{s}(b)={\bf x}_{t}(b))$. First, we note that when the pairs $(i,j)$ and $(s,t)$ are disjoint, the events ${\bf x}_{i}(a)={\bf x}_{j}(a)$ and ${\bf x}_{s}(b)={\bf x}_{t}(b)$ are independent:

$\displaystyle\Pr({\bf x}_{i}(a)={\bf x}_{j}(a),{\bf x}_{s}(b)={\bf x}_{t}(b))$

$\displaystyle=n^{-2\beta}.$

The events are also independent when $i=s$ but $j\neq t$, or $j=t$ but $i\neq s$,

$\displaystyle\Pr({\bf x}_{i}(a)={\bf x}_{j}(a),{\bf x}_{i}(b)={\bf x}_{t}(b))=n^{-2\beta},$

and when $i=s$, $j=t$, and $|a-b|\geq k$. However, when $i=s$, $j=t$, and $|a-b|<k$,

$\displaystyle\Pr({\bf x}_{i}(a)={\bf x}_{j}(a),{\bf x}_{i}(b)={\bf x}_{j}(b))=2^{-(k+|a-b|)}\leq n^{-\beta}$

as there are $|a-b|$ symbols of overlap between ${\bf x}_{i}(a)$ and ${\bf x}_{j}(b)$. Therefore,

	$\displaystyle\mathbb{E}[V^{2}]$	$\displaystyle=\mathbb{E}\left[\sum_{i=1}^{m-1}\sum_{j=i+1}^{m}\sum_{a=1}^{n-k+1}Z_{i,j,a}\sum_{s=1}^{m-1}\sum_{t=s+1}^{m}\sum_{b=1}^{n-k+1}Z_{s,t,b}\right]$
		$\displaystyle\leq m^{4}n^{2-2\beta}+2m^{2}n^{1-\beta}k.$		(10)

Combining $\mathbb{E}[V]$ and (10), we have that

	$\displaystyle\Pr({\mathcal{E}})$	$\displaystyle\geq\frac{\binom{m}{2}^{2}(n-k+1)^{2}n^{-2\beta}}{m^{4}n^{2-2\beta}+2m^{2}n^{1-\beta}k}$		(11)
		$\displaystyle\sim\frac{m^{4}n^{2}n^{-2\beta}}{4m^{4}n^{2-2\beta}+8m^{2}n^{1-\beta}k}$		(12)
		$\displaystyle=\frac{1}{4+\frac{8\beta\log n}{n^{2\alpha+1-\beta}}}$		(13)

where (12) is asymptotically equivalent as $n\to\infty$. Finally, we notice that if $\beta<2\alpha+1$, (13) converges to $\frac{1}{4}$, implying that $\Pr({\mathcal{E}})$ is bounded away from zero and $(\alpha,\beta)\notin{\mathcal{F}}$.

Next, we show that $\Pr({\mathcal{E}})\not\to 0$ when $\beta<\alpha+\frac{3}{2}$. For this regime, we will define the event ${\mathcal{H}}$ that there are two strings ${\bf x}_{i}$, ${\bf x}_{j}$ and indices $a$, $b$, $c$ such that

$${\bf x}_{i}(a)={\bf x}_{j}(c),{\bf x}_{i}(b)={\bf x}_{j}(c+b-a).$$

This is illustrated in Figure 2(a). Essentially, we have two repeated $k$-mers across sequences ${\bf x}_{i}$ and ${\bf x}_{j}$, and the gap between them is the same in both sequences. Given ${\mathcal{H}}$, an incorrect reconstruction exists because we can swap the middle parts of ${\bf x}_{i}$ and ${\bf x}_{j}$, creating

	$\displaystyle\tilde{\bf x}_{i}={\bf x}_{i}[1:a-1]{\bf x}_{j}[c:c+b-a-1]{\bf x}_{i}[b:n]$
	$\displaystyle\tilde{\bf x}_{j}={\bf x}_{j}[1:c-1]{\bf x}_{i}[a:b-1]{\bf x}_{j}[c+b-a:n],$

which can be verified to be two length-$n$ sequences. We will show that $\Pr({\mathcal{H}})\not\to 0$ as $n\to\infty$.

Define $W_{i,j,\mathbf{a}}=\mathbbm{1}_{{\bf x}_{i}(a_{1})={\bf x}_{j}(a_{3}),{\bf x}_{i}(a_{2})={\bf x}_{j}(a_{3}+a_{2}-a_{1})}$ for $\mathbf{a}=(a_{1},a_{2},a_{3})$ and let $U=\sum_{i,j}\sum_{\mathbf{a}}W_{i,j,\mathbf{a}}$. We can again use the Paley-Zygmund inequality to find a lower bound on $\Pr(U>0)$. The second moment $\mathbb{E}[U^{2}]$ can be computed similarly to $\mathbb{E}[V^{2}]$, considering separately the terms of the form $\mathbb{E}[W_{ij\mathbf{a}}W_{ij\mathbf{b}}]$ where the $k$-mers indicated by $\mathbf{a}$ overlap with those indicated by $\mathbf{b}$:

	$\displaystyle\mathbb{E}[U^{2}]$	$\displaystyle=\mathbb{E}\left[\sum_{i=1}^{m-1}\sum_{j=i+1}^{m}\sum_{\mathbf{a}}W_{i,j,\mathbf{a}}\sum_{s=1}^{m-1}\sum_{t=s+1}^{m}\sum_{\mathbf{b}}W_{s,t,\mathbf{b}}\right]$
		$\displaystyle\leq(m^{2}n^{3}n^{-2\beta})^{2}+m^{2}n^{3}(2k)^{3}n^{-2\beta}$		(14)

where the first term in (14) covers all $k$-mers that are independent of each other and the second term covers the case where $i=s$, $j=t$, and $a_{i}$ overlaps with $b_{i}$ for $i=1,2,3$.

To compute the first moment, $\mathbb{E}[U]$, we must account for which assignments of $\mathbf{a}$ are possible. More precisely, we need to count the vectors $(a_{1},a_{2},a_{3})\in[1:n^{\prime}]^{3}$ such that $a_{1}<a_{2}$ and $a_{2}+a_{3}-a_{1}\in[1:n^{\prime}]$, where $n^{\prime}=n-k+1$. This is cumbersome, so we consider a lower bound by considering $a_{1}\in[1:n^{\prime}/4]$, $a_{2}\in[n^{\prime}/4+1:n^{\prime}/2]$, and $a_{3}\in[n^{\prime}/4+1:n^{\prime}/2]$. Notice that for any of the $(n^{\prime}/4)^{3}$ choices, $a_{2}+a_{3}-a_{1}\in[1:n^{\prime}]$. We obtain the (rather loose) lower bound

$\displaystyle\mathbb{E}[U]\geq\binom{m}{2}\left(\frac{n^{\prime}}{4}\right)^{3}n^{-2\beta}.$

(15)

Therefore, we can lower bound the error probability as

	$\displaystyle\Pr({\mathcal{E}})$	$\displaystyle\geq\frac{\mathbb{E}[U]^{2}}{\mathbb{E}[U^{2}]}\geq\frac{\binom{m}{2}^{2}(n^{\prime}/4)^{6}n^{-4\beta}}{m^{4}n^{6}n^{-6\beta}+m^{2}n^{3}(2k)^{3}n^{-2\beta}}$
		$\displaystyle\sim\frac{4^{-7}m^{4}n^{6-4\beta}}{m^{4}n^{6-4\beta}+8m^{2}k^{3}n^{3-2\beta}}$
		$\displaystyle=\frac{4^{-7}}{1+\frac{8k^{3}}{n^{2^{\alpha+3-2\beta}}},}$

which tends to $4^{-7}$ as $n\to\infty$, as long as $\beta<\alpha+\tfrac{3}{2}$. We conclude that for any $(\alpha,\beta)$ such that $\beta<\alpha+\tfrac{3}{2}$ or $\beta<2\alpha+1$, $\Pr({\mathcal{E}})\not\to 0$ as $n\to\infty$. This concludes the proof of the outer bound in Theorem 1.

In this paper, we have characterized the feasibility of a multiple sequence reconstruction problem using $k$-mers for a large region of $(\alpha,\beta)$ parameters. We identified a region where the $k$-mers of each string are disjoint with high probability and correct recovery of the source sequences is trivial, as well as a region where, while repeats exist, the only possible reconstruction of the source sequences is the correct one. In the latter repeat-abundant region, careful examination of the de Bruijn graph $G({\mathcal{Y}})$ reveals the unique set of paths covering the graph. Finally, we have determined a region where incorrect reconstructions exist with positive probability.

The above findings have characterized nearly all $(\alpha,\beta)$ pairs for $\alpha>0$, $\beta>0$, except for a small region between $\beta=\max(2\alpha+1,\alpha+\frac{3}{2})$ and $\beta=\alpha+2$. We believe that part of the remaining region may be where the solution is not unique, and close analysis of additional difference de Bruijn graph topologies for a larger number of strings and shared $k$-mers might result in a larger red region in Figure 1.

This work invites future study in several directions. One of these is generalization to other source distributions; some results here may be immediately trivially extended to a general i.i.d. source. Additionally, motivated by real-world DNA sequencing technologies, future work may analyze source reconstruction from noisy $k$-mers or a random number of copies of each $k$-mer.

The work of I.S. was supported in part by the National Science Foundation under CCF grants 2007597 and 2046991. The authors would like to thank the anonymous reviewers of this manuscript for their detailed feedback and suggestions for improvement.