Capacity Results for the Noisy Shuffling Channel

TL;DR

This paper analyzes the capacity of a noisy shuffling channel inspired by DNA storage, providing exact capacity results for the BSC-shuffling model and demonstrating the optimality of a simple index-based coding scheme.

Contribution

It characterizes the capacity of the noisy shuffling channel and proves the optimality of a straightforward coding scheme, advancing understanding of DNA storage channel limits.

Findings

Exact capacity characterization for the BSC-shuffling channel

Optimality of a simple index-based coding scheme

Relevance to DNA storage systems

Abstract

Motivated by DNA-based storage, we study the noisy shuffling channel, which can be seen as the concatenation of a standard noisy channel (such as the BSC) and a shuffling channel, which breaks the data block into small pieces and shuffles them. This channel models a DNA storage system, by capturing two of its key aspects: (1) the data is written onto many short DNA molecules that are stored in an unordered way and (2) the molecules are corrupted by noise at synthesis, sequencing, and during storage. For the BSC-shuffling channel we characterize the capacity exactly (for a large set of parameters), and show that a simple index-based coding scheme is optimal.

Equations

$$(1-e^{-c})(1-1/\beta ),$$

$$C=1-H(p)-1/\beta .$$

$$X^{ML}=\left[X_1^L,X_2^L,\dots ,X_M^L\right]$$

$$Y^{ML}=\left[Y_1^L,Y_2^L,\dots ,Y_M^L\right]$$

$$Y_k^L$$

$$\frac{M\left(L{R}_{\text{BSC}}-\log M\right)}{ML}=R_{\text{BSC}}-1/\beta .$$

$$R_{\text{index}}=1-H(p)-1/\beta ,$$

$$C\le \min \left[1-H(p),1-1/\beta \right],$$

$$C=1-H(p)-1/\beta .$$

$$X^{ML}=\left[X_1^L,X_2^L,\dots ,X_M^L\right]$$

$$Y^{ML}=\left[Y_1^L,Y_2^L,\dots ,Y_M^L\right]$$

$$MLR$$

$$ML(R-{ϵ}_M)$$

$$=H\left(Y^{ML}\right)-H\left(S^M,Z^{ML},Y^{ML}|X^{ML}\right)$$

$$+H\left(S^M,Z^{ML}|X^{ML},Y^{ML}\right)$$

$$=H\left(Y^{ML}\right)-H\left(S^M,Z^{ML},Y^{ML}|X^{ML}\right)$$

$$+H\left(S^M|X^{ML},Y^{ML}\right).$$

$$Z_{S(k)}^L=Y_k^L\oplus X_{S(k)}^L$$

$$H\left(S^M,Z^{ML},Y^{ML}|X^{ML}\right)$$

$$=H\left(S^M\right)+H\left(Z^{ML}\right)+H\left(Y^{ML}|X^{ML},S^M,Z^{ML}\right)$$

$$=\log M!+MLH(p)+0$$

$$=M\log M+MLH(p)+o(ML).$$

$$1-H(2p)-\frac{2}{\beta }>0$$

$$H\left(Y^{ML}\right)+H\left(S^M|X^{ML},Y^{ML}\right)\le ML+o(ML).$$

$$ML(R-{ϵ}_M)\le ML-MLH(p)-M\log M+o(ML).$$

$$H\left(Y^{ML}\right)+H\left(S^M|X^{ML},Y^{ML}\right)\le ML+M\log M.$$

$$H\left(Y^{ML}\right)\le M+L\mathbf{E}|T|$$

$$+(M-\mathbf{E}|T|)\left(\log \mathbf{E}|T|+LH(\alpha )\right),$$

$$H\left(S^M|X^{ML},Y^{ML}\right)\le M\log M$$

$$-\mathbf{E}|T|\log \mathbf{E}|T|+o(ML).$$

$$H\left(Y^{ML}\right)+H\left(S^M|X^{ML},Y^{ML}\right)$$

$$\le L\mathbf{E}|T|+(M-\mathbf{E}|T|)\left(\log \mathbf{E}|T|+LH(\alpha )\right)$$

$$+M\log M-\mathbf{E}|T|\log \mathbf{E}|T|+o(ML)$$

$$=\mathbf{E}|T|L(1-H(\alpha ))+M\log \mathbf{E}|T|-2\mathbf{E}|T|\log \mathbf{E}|T|$$

$$+MLH(\alpha )+M\log M+o(ML).$$

$$f(x)\triangleq x\beta (1-H(\alpha ))\log M+M\log x-2x\log x.$$

$$f^{\prime }(x)$$

Full text

I Introduction

Due to its longevity and enormous information density, and thanks to rapid advances in technologies for writing (synthesis) and reading (sequencing), DNA is on track to become an attractive medium for archival data storage. Computer scientists and engineers have dreamed of harnessing DNA’s storage capabilities already in the 60s [1, 2], and in recent years this idea has developed into an active field of research: In 2012 and 2013 groups lead by Church [3] and Goldman [4] independently stored about a megabyte of data in DNA. In 2015 Grass et al. [5] demonstrated that millenia long storage times are possible by protecting the data both physically and information-theoretically, and designed a robust DNA data storage scheme using modern error correcting codes. Later, in the same year, Yazdi et al [6] showed how to selectively access files, and in 2017, Erlich and Zielinski [7] demonstrated that practical DNA storage can achieve very high information densities. In 2018, Organick et al. [8] scaled up these techniques and stored about 200 megabytes of data.

DNA is a long molecule made up of four nucleotides (Adenine, Cytosine, Guanine, and Thymine) and, for storage purposes, can be viewed as a string over a four-letter alphabet. However, due to technological constraints, it is difficult and inefficient to synthesize long strands of DNA. Thus, data is stored on short DNA molecules which are preserved in a pool of DNA and cannot be spatially ordered. Indeed, all the aforementioned works did store information on molecules of no longer than a few hundred nucleotides.

Given these constraints, the simplest mathematical model for a DNA storage system is a shuffling channel: data is encoded into $M$ DNA molecules (i.e., strings), each of length $L$, and the output of the channel consists of a random shuffle of the $M$ sequences.

In our previous work [9], we studied the capacity of such a shuffling channel with one additional feature: random sampling of the shuffled sequences. In a DNA storage system, the information is accessed via shotgun sequencing, which in essence corresponds to randomly sampling and reading molecules from the DNA pool. In order to formally study the capacity of this system, we considered an asymptotic regime in which $M$ is the number of sequences/molecules stored, $N=cM$ sequences are sampled, the length of each sequence is small (specifically $L=\beta\log M$ for a fixed $\beta>0$, and throughout $\log$ is with respect to base $2$), and $M\to\infty$. For this setting, and assuming an alphabet of size two for simplicity, the storage capacity, defined as the maximum number of bits that can be reliably stored per symbol (the total number of symbols is $ML$) was shown to be

[TABLE]

provided that $\beta>1$. If $\beta<1$, no positive rate is achievable.

Surprisingly, the capacity of this shuffle-and-sampling channel is achieved by a simple scheme where each sequence is given a unique index (or header), and an erasure outer code is used to protect against missing molecules in the sampling step. The factor $1-e^{-c}$ can be understood as the loss due to unseen molecules, and the factor $1-1/\beta$ corresponds to the loss due to shuffling or the lack of ordering in the molecules.

In practice, however, the DNA sequences are additionally corrupted by substitutions, deletions, and insertion errors (which occur both at synthesis, sequencing, and during the storage period). For that reason, several recent works have studied the design of error-correcting codes tailored to the error profiles of DNA storage systems [10, 6, 11, 12, 13]. However, generalizing the capacity expression in (1) to a setting in which, in addition to random shuffling and random sampling, the molecules are corrupted by noise is quite nontrivial.

In this work we make progress on this front by studying the capacity of a noisy shuffling channel. We drop the random sampling constraint considered in [9], and focus on a setting in which the $M$ sequences are first corrupted by noise and then randomly shuffled. Since in all recently proposed systems [3, 4, 5, 6, 7, 8], substitution errors are more frequent than insertion and deletion errors [14], we focus on substitution noise.

Contributions: We study the capacity of the noisy shuffling channel illustrated in Figure I, where the input/ouput sequences are binary, and the noisy channel is a Binary Symmetric Channel (BSC) with crossover probability $p$. Our main result establishes that, for a large set of the parameters $\beta$ and $p$, the capacity of the BSC shuffling channel is

[TABLE]

Since the capacity of a BSC channel is $C_{\text{BSC}}=1-H(p)$, this results implies that the shuffling reduces the BSC capacity by $1/\beta$ bits.

The capacity-achieving scheme consists of encoding each length-$L$ string as a length-$L$ codeword from a BSC capacity-achiving code, and using $\log M$ out of its $LC_{\text{BSC}}$ information bits (prior to encoding) to store an index. Hence, at least asymptotically, there is no rate gain in coding across the $M$ strings, and index-based approaches, which were used in several practical implementations [5, 6, 7, 8] are optimal.

**Outline: ** The paper is organized as follows. After formalizing the problem setting in Section II, we state our main result and describe the achievability in Section III. Section IV is dedicated to the converse. We conclude by discussing extensions and future directions in Section V.

**Related literature: ** Motivated by DNA-based storage, a few recent works have considered the problem of coding across an unordered set of strings [15, 16, 17, 18]. The setting studied in all these works bears similarities with the one in this paper, but they focus on providing explicit code constructions, as opposed to characterizing the channel capacity, as we do here.

II Problem Setting

The BSC shuffling channel, illustrated in Figure I, can be seen as the concatenation of a BSC with crossover probability $p$ and a shuffling channel, which randomly reorders the binary strings. Formally, we will view the input to the channel as a length-$ML$ binary string

[TABLE]

(i.e., $M$ strings of length $L$ concatenated to form a single string of length $ML$). Similarly,

[TABLE]

is the corresponding output. We let $Z^{ML}=\left[Z_{1}^{L},\ldots,Z_{M}^{L}\right]$ be the random binary error pattern created by the BSC on $X^{ML}$ and $S^{M}\in\{1,\ldots,M\}^{M}$ be a vector encoding the random shuffle induced by the channel, so that

[TABLE]

for $k=1,\ldots,M$, where $\oplus$ indicates elementwise modulo $2$ addition. Similar to our previous work [9], we let $L=\beta\log M$, and consider the asymptotic regime $M\to\infty$. As shown in [9], this is the proper regime for positive data rates to be achievable.

A rate $R$ is said to be achievable if we have a sequence of codes, each with $2^{MLR}$ codewords, and whose error probability goes to zero as $M\to\infty$. The shuffling channel capacity $C$ is the supremum over all achievable rates.

III Main Result

Our main result establishes that, provided that $\beta$ is large enough, treating each length-$L$ sequence as the input to a separate BSC and encoding a unique index into each sequence is capacity optimal. More precisely, consider code for a BSC with codewords of length $L$ and rate $R_{\text{BSC}}=1-H(p)-\epsilon$, for some small $\epsilon$. Using this code, we can store $LR_{\text{BSC}}$ information bits in each length-$L$ string. We utilize the first $\log M$ of those bits to encode a distinct index for each of the $M$ strings. Hence, in each string we can effectively store $LR_{\text{BSC}}-\log M$ data bits, and the rate of this scheme is

[TABLE]

Since $\epsilon>0$ can be chosen arbitrarily small, this scheme can have a rate arbitrarily close to

[TABLE]

if $1-H(p)-1/\beta>0$.

Strictly speaking, the simple index-based scheme described above needs to be modified using an outer code to guarantee that the rate in (3) is actually achieved. That is because there is a vanishing, but positive, probability that a given string is decoded in error, and since $M\to\infty$, the probability that at least one of the strings is decoded in error may not go to zero. When an error occurs, the index is also read in error, which may cause a “collision” with a properly decoded string. Nevertheless, since the error probability for each given string vanishes as $M\to\infty$, this can be fixed with an outer code across the data symbols of all the $M$ strings, and it is straightforward to see that one can achieve arbitrarily close to the rate in (3). Hence, (4) is a lower bound to the capacity $C$.

On the other hand, a simple upper bound is given by

[TABLE]

since both the capacity of the BSC, $C_{\text{BSC}}=1-H(p)$, and the capacity of the shuffling channel, $1-1/\beta$ [9], must be upper bounds to $C$.

Our main result improves on the upper bound in (5), and establishes that for parameters $(p,\beta)$ in a certain regime, the lower bound in equation (4) is actually the capacity.

Theorem 1**.**

If $p\leq 0.1$ and $\beta\geq 6.4$,

[TABLE]

*Moreover, if $\beta\leq 1$, $C=0$. *

As we will see, the converse argument presented in the next section and the capacity expression in (6) hold for a set of parameters $(p,\beta)$ larger than those specified in Theorem 1. In fact, for all $(p,\beta)$ in the blue region of Figure III, the capacity is given by (6). Hence, for $p<0.01$ for instance, we need $\beta\geq 2.35$ for (6) to hold.

IV Converse

Consider a sequence of codes for the noisy shuffling channel with rate $R$ and vanishing error probability. Let

[TABLE]

be the input to the channel when we choose one of the $2^{MLR}$ codewords from one such code uniformly at random, and

[TABLE]

the corresponding output. We first observe that

[TABLE]

where $\epsilon_{M}\to 0$ as $M\to\infty$ by Fano’s inequality. Then,

[TABLE]

The last equality follows by noticing that, given $(S^{M},X^{ML},Y^{ML})$, one can compute

[TABLE]

for $k=1,\ldots,M$, and thus $H\left(Z^{ML}|X^{ML},Y^{ML},S^{M}\right)=0$. The second term in (7) can be computed as

[TABLE]

In order to finish the converse, we need to jointly bound the first and third terms in equation (7). This is the most challenging part of the proof and is summarized in the following lemma:

Lemma 1**.**

If $\beta$ and $p$ satisfy

[TABLE]

then it holds that

[TABLE]

The parameter regime $(p,\beta)$ for which the inequality (9) holds is the regime in which our capacity expression holds. The regime is depicted in blue in Figure III.

Combining (7), (8) and Lemma 1, we obtain

[TABLE]

Dividing by $ML$ and letting $M\to\infty$ yields the converse.

IV-A Intuition for Lemma 1

If we trivially bound each entropy term separately, we obtain

[TABLE]

However, intuitively, the bound $H\left(S^{M}|X^{ML},Y^{ML}\right)\leq M\log M$ is too loose. Given $X^{ML}=x^{ML}$ and $Y^{ML}=y^{ML}$, one could estimate $S^{M}$ by finding for each $y_{i}^{L}$ the $x_{j}^{L}$ that is closest to it and guessing that $S(i)=j$. This will be a good guess if there is no other $x_{k}^{L}$ close to $x_{j}^{L}$. Hence, we can picture two opposite scenarios, illustrated in Figure IV-A.

In the first case, the strings $x_{1}^{L},\ldots,x_{M}^{L}$ are all distant from each other (in the Hamming sense). Hence, the maximum likelihood estimate of $S^{M}$ given $X^{ML}=x^{ML}$ and $Y^{ML}=y^{ML}$ will be “close” to the truth and we would expect $H\left(S^{M}|X^{ML}=x^{ML},Y^{ML}=y^{ML}\right)$ to be small. In the second case, illustrated in Fig. IV-A(b), several of $x_{1}^{L},...,x_{M}^{L}$ are close to each other. So we have less information about $S^{M}$, and $H\left(S^{M}|X^{ML}=x^{ML},Y^{ML}=y^{ML}\right)$ may be large.

The term $H\left(Y^{ML}\right)$ is maximized by making the $\left\{X_{i}^{L}\right\}$ independent and taking values uniformly in $\{0,1\}^{L}$. Hence, in order for $H\left(Y^{ML}\right)$ to be large, we would expect the scenario in Fig. IV-A(a) to be common and the scenario in Fig. IV-A(b) to be rare. This creates a tension between the terms $H\left(Y^{ML}\right)$ and $H\left(S^{M}|X^{ML},Y^{ML}\right)$, which we exploit to prove Lemma 1.

IV-B Proof of Lemma 1

Let $\{Y^{L}_{1},\ldots,Y^{L}_{M}\}$ be the set of $M$ strings observed at the output of the channel. In order to capture whether we are in the scenario of Figure IV-A(a) or (b), we let $T$ be the largest subset of $[1:M]\triangleq\{1,\ldots,M\}$ so that, for any $i,j\in T$, $d_{H}\left(Y_{i}^{L},Y_{j}^{L}\right)\geq\alpha L,$ where $d_{H}$ is the Hamming distance and $\alpha>2p$. We assume that in case of ties, an arbitrary tie-breaking rule is used to define $T$ (the actual choice will not be relevant for the proof).

We will prove that, given the conditions in Lemma 1, the following two bounds involving $\mathbf{E}|T|$ hold:

[TABLE]

For large $\mathbf{E}|T|$, we are intuitively more often in the scenario in Figure IV-A(a), while Figure IV-A(b) corresponds to the case where $\mathbf{E}|T|$ is small. The bounds above capture the tension between $H\left(Y^{ML}\right)$ and $H\left(S^{M}|X^{ML},Y^{ML}\right)$ because (B2) is decreasing in $\mathbf{E}|T|$, while (B1) is increasing in $\mathbf{E}|T|$ (as long as $\beta(1-H(\alpha))\geq 1$). Combining (B1) and (B2),

[TABLE]

Replacing $\mathbf{E}|T|$ with $x$ and ignoring the terms in this upper bound that do not involve $x$, we have the expression

[TABLE]

Let $\gamma=\beta(1-H(\alpha))$. For $x>0$, we have

[TABLE]

Hence $f^{\prime}(x)>0$ if

[TABLE]

We see that, as long as $\gamma>2$, the right-hand side of (11) is greater than $M$ for $M$ large enough. This means that $f(x)$ is increasing for $1\leq x\leq M$, and must attain its maximum at $f(M)$. Notice that $x=\mathbf{E}|T|\leq M$. Therefore, (10) can be upper-bounded by setting $x=\mathbf{E}|T|=M$, which yields

[TABLE]

Notice that this holds as long as

[TABLE]

for some $\alpha>2p$. From the continuity of $H(\cdot)$, such $\alpha$ can be found if (9) is satisfied, proving the lemma. Finally, we prove (B1) and (B2).

*Proof of (B1): * Since $T$ is a deterministic function of $Y^{ML}$ and can take at most $2^{M}$ values,

[TABLE]

Next we notice that, for a given $t$, we can write

[TABLE]

The first term in (13) is trivially bounded as

[TABLE]

Each of the remaining length-$L$ strings $Y^{L}_{i}$ with $i\notin t$ must be within a distance $\alpha L$ from one of the strings in $\{Y^{L}_{i}:i\in t\}$, from the definition of $T$. Hence, conditioned on $\{Y^{L}_{i}:i\in t\}$, each of them can only take at most $|t||B(\alpha L)|$ values, where $B(\alpha L)$ is a Hamming ball of radius $\alpha L$. Since $|B(\alpha L)|\leq 2^{LH(\alpha)}$ for $\alpha<0.5$, we bound the second term in (13) as

[TABLE]

Using these bounds back in (12), we obtain

[TABLE]

where we used the fact that $(M-x)\log x$ is a concave function of $x$ and Jensen’s inequality. ∎

*Proof of (B2): * Since $T$ is a deterministic function of $Y^{ML}$,

[TABLE]

Next we notice that the probability that $\delta L$ or more errors occur in a single length-$L$ string, for $\delta>p$, is at most $2^{-LD(\delta||p)}$ by the Chernoff bound (where $D(\cdot||\cdot)$ is the binary KL divergence). If we let ${\mathcal{E}}_{i}$ be the event that $d_{H}\left(X_{S(i)}^{L},Y_{i}^{L}\right)\geq\delta L$, then we have

[TABLE]

The conditional entropy term in (15) is upper bounded by

[TABLE]

The final step is to bound the conditional entropy term in (16), for the case where $i\in t$. Set $\delta=\alpha/2$. Conditioned on $\bar{\mathcal{E}}_{i}$, $d_{H}\left(X_{S(i)}^{L},Y_{i}^{L}\right)<\alpha L/2$. Moreover, conditioned on $T=t$, for any $j\in t-\{i\}$, $d_{H}\left(Y_{i}^{L},Y_{j}^{L}\right)\geq\alpha L$. Define the set

[TABLE]

We claim that, if $i\in t$, $S(i)$ must be in $A_{i}$. To see this notice that, for any $k\in t$, $k\neq i$, we have

[TABLE]

implying that $d_{H}\left(X_{S(i)}^{L},Y_{k}^{L}\right)>\alpha L/2\geq d_{H}\left(X_{S(i)}^{L},Y_{i}^{L}\right)$. Therefore, $S(i)$ for each output string $Y^{L}_{i}$ with $i\in t$, can take at most $|A_{i}|$ values. Hence we have

[TABLE]

where the last inequality follows because $\sum_{i\in t}|A_{i}|=M$, and the sum is maximized by $|A_{i}|=M/|t|$. Combining (15), (16), and (17), we obtain

[TABLE]

where we used Jensen’s inequality in the last step. Since $M^{-\beta D(\delta||p)}\to 0$ as $M\to\infty$, $M^{-\beta D(\delta||p)}M\log M=o(ML)$, concluding the proof. ∎

V Discussion and Extensions

While the parameter regime of $(p,\beta)$ in the blue region of Figure III is practically the most relevant one, an interesting question for future work is whether expression (6) is still the capacity of the BSC-shuffling channel if $\beta$ and $p$ do not satisfy (9) (i.e., the gray region in Figure III). Notice that this is a high-noise, short-block regime, and it is reasonable to imagine that coding across the different sequences can be helpful and an index-based approach might not be optimal.

V-A General symmetric channels

If we view the capacity expression in (6) as $C_{\text{BSC}}-1/\beta$, it is natural to wonder whether for a different noisy channel with capacity $C_{\text{noisy}}$, the correspoinding noisy shuffling channel would have capacity $C_{\text{noisy}}-1/\beta$. As it turns out, the converse proof in Section IV can be extended to the class of symmetric discrete memoryless channels (as described in [19, ch. 7.2]).

Consider a noisy shuffling channel formed by a symmetric discrete memoryless channel (SDMC) with output alphabet $\mathcal{Y}$, followed by a shuffling channel. We focus on the same asymptotic regime from Section II. Then we have:

Theorem 2**.**

If $\beta$ is large enough, the capacity of the SDMC-shuffling channel is given by

[TABLE]

Moreover, if $\beta\leq\log|\mathcal{Y}|$, $C=0$.

For symmetric channels, capacity is achieved by making the distribution of the output uniform, which allows an analogous result to Lemma 1 to be obtained. Notice that how large $\beta$ needs to be depends on the specific channel transition matrix.

VI Concluding Remarks

In this paper we took steps towards the understanding of the fundamental limits of DNA-based storage systems. We proposed a simple model capturing the fact that molecules are stored in an unordered fashion, are short, and are corrupted by individual base errors. Our results show that a simple index-based coding scheme is asymptotically optimal.

While the model captures (moderate) substitution errors which are the prevalent error source on a nucleotide level of current DNA storage systems, the current generation of systems relies on low-error synthesis and sequencing technologies that are relatively expensive and limited in speed. A key idea towards developing the next-generation of DNA storage systems is to employ high-error, but cheaper and faster synthesis and sequencing technologies such as light-directed maskless synthesis of DNA and nanopore sequencing. Such systems induce a significant amount of insertion and deletion errors. Thus, and important area of further investigation is to understand the capacity of channels which introduce deletions and insertions as well.

In addition, in practice, one can introduce physical redundancy (i.e., duplicate molecules) in storage and redundancy in the DNA sequencing. Thus, an interesting question is to understand how to efficiently code for channels where we have several unordered, independently perturbed copies of each sequence.