On Coding for an Abstracted Nanopore Channel for DNA Storage

In the emerging field of DNA storage, data is encoded as DNA sequences and stored; the data can be read back by sequencing the stored DNA. This technology promises high storage density and stability, as well as efficient duplication of data and random access using PCR-based technologies; we refer the reader to [4] and the references therein for an excellent overview.

Both the synthesis and sequencing processes are noisy, and as a result the data must be encoded before the synthesis stage to ensure accurate data recovery. Prior work has studied methods for encoding data in order to protect it against (aspects of) the noise introduced by these processes, for example [10, 3, 7, 22, 17, 14, 15, 5].

In this work, we focus on one particular stage of this noisy process, the nanopore sequencer. Nanopore sequencing—and in particular the MinION sequencer developed by Oxford Nanopore Technologies [12]—is an emerging sequencing technology. While initial works in DNA storage used Illumina sequencing, nanopore sequencing has been attracting interest due to its portability, low cost, and ability to support significantly longer reads than Illumina.

At a high level, the nanopore sequencer works as follows. A single strand of DNA is passed through a pore, leading to variations in a current readout. The pore can hold $k$ nucleotides (for our purposes, a nucleotide is just a value in $\{A,C,G,T\}$) at a time; in practice $k$ is about six. The value of the current readout depends on which nucleotides are in the pore. For example, if the strand of DNA is AGCTAAT, and the pore sees the sub-strand AGCT, it will output one current reading. As the strand is passed through the pore, the contents of the pore will shift, say from AGCT to GCTA, then to CTAA, then to TAAT, and so on. This will result in a change in the current reading, according to some function $f$ that maps $k$-mers to current levels. The process is depicted in Figure 1. Given the current readout, the goal is to recover the original DNA sequence.

This channel is difficult to analyze, for several reasons. First, the output at any given time depends on $k>1$ bases, and so there is inter-symbol interference. Second, there may be collisions in the output: two different pore contents may lead to similar current readouts. Third, the current readout can be noisy. Fourth, the amount of time that each $k$-mer spends in the pore can vary, and sometimes never occur at all, leading to synchronization errors in the output.

Due to this complexity, typical basecallers (that is, methods for recovering the original sequence from the current readouts) rely on machine learning techniques [21, 9, 20]. This is effective in practice, but difficult to get a theoretical handle on. While there are several practical approaches to error correction for nanopore sequencing in DNA storage [22, 17, 14, 15, 5], the theoretical limits of this channel are not well understood. The work [16], discussed more below, proposed a probabilistic model of the nanopore sequencer and developed bounds on the capacity of the resulting channel. However, this problem is quite difficult; in particular, taking into account the possibility of synchronization errors places this problem as more difficult than determining the capacity of the deletion channel, a notorious open problem.

In this work, we take a different approach to developing a theoretical understanding of the nanopore sequencer. We develop a highly abstracted, deterministic model, which is meant to highlight the first two sources of noise mentioned above: the inter-symbol interference, and the prospect of collisions when multiple $k$-mers lead to similar current outputs. We develop methods and theory to understand the capacity of our abstracted channel, and we propose efficient coding schemes and algorithms.

Our goal in this work is to open up a research direction towards a theoretical understanding of the nanopore sequencer. We fully acknowledge that our work is not yet practical; in particular our abstracted channel does not include noise in the current readings or synchronization errors that can arise from variable pore dwelling times. It is our hope that a solid understanding of our abstracted channel can then be combined with (much more well-studied) theories of coding for substitution and synchronization errors, in order to make progress on a more realistic channel model.

Contributions. Our contributions are as follows.

1. 1.

   We propose a novel abstraction of the nanopore sequencer, which highlights the inter-symbol interference and the prospect of collisions. This abstraction is simple enough that it is tractable, yet complex enough that it (a) captures some fundamental properties of the nanopore sequencer, and (b) already gives rise to extremely interesting problems from a theoretical perspective. We hope that our abstraction will lead to future work in this area.

2. 2.

   In Section 3, we develop an algorithm to determine the capacity of this channel. The algorithm is inefficient as pore size $k$ grows, but we can use it to determine the capacity for small $k$.

3. 3.

   In Section 4, we develop simple bounds on the “best” and “worst” capacity for an arbitrary pore size $k$ (where “best” and “worst” refers to the choice of the map from $k$-mers to current readouts). We use our algorithms mentioned above to compare these bounds to the exact values for $k=2$. These preliminary results suggest there may be some intriguing “bumpiness” in the worst-case capacity as the number of distinct current levels increases, but that the average-case is closer to best- than worst-case.

4. 4.

   In Section 5, we develop efficient coding schemes for our abstracted channel. In particular, we develop a scheme that achieves rate $C_{f}-\epsilon$, where $C_{f}$ is the capacity, that requires preprocessing time $O\left(\frac{1}{\epsilon}\cdot 4^{1/\epsilon}\right)$, and has encoding/decoding time $O(n)$, where the $O(\cdot)$ notation hides a constant that depends on $f$. We also present an efficient coding scheme that improves over the “naive” one with high probability, where the probability is over a random map from $k$-mers to current readouts.

In this section we formalize our model. As mentioned above, our goal is to focus on the challenges of (1) inter-symbol interference, and (2) the possibility of different $k$-mers producing similar current readouts. With that in mind, we propose a very simple model for the nanopore sequencer. As input, we take an encoded string $s_{0}s_{1}\cdots s_{n-1}\in\{A,C,G,T\}^{n}$. This is transformed into a sequence of $k$-mers according to a sliding window, to obtain $(s_{0}\cdots s_{k-1}),(s_{1}\cdots s_{k}),(s_{2}\cdots s_{k+1}),\ldots,(s_{n-k}\cdots s_{n-1})$. Finally, each of these $k$-mers is mapped to one of $b$ distinct current levels, according to a mapping $f:\{A,C,G,T\}^{k}\to\{0,1,\ldots,b-1\}$. This mapping $f$ defines the channel.

Given a mapping $f$, we are interested in the capacity $C_{f}$ (in bits-per-base), which we define as follows.

Observe that if $\mathcal{S}$ is a collection of strings $s\in\{A,C,G,T\}^{n}$ so that $f^{*}(s)$ are all distinct for $s\in\mathcal{S}$, then by assigning a different message to each $s\in\mathcal{S}$, we can communicate perfectly (if not necessarily efficiently) across the abstract nanopore channel.

We will focus on balanced mappings $f$. These are mappings so that $|f^{-1}(i)|$ is the same for all $i\in\{0,1,\ldots,b-1\}$.

Our first contribution is an algorithm (Algorithm 1 below) that computes $C_{f}$, given $f$. The basic idea is to consider a finite automaton on the alphabet $\{0,1,\dots,b-1\}$ that accepts exactly those current readouts that can be generated by some DNA strand; then we use the transfer matrix method [2, 8] for counting accepting paths in that finite automaton.

Formally, an Nondeterministic Finite Automaton (NFA) is a tuple $(Q,\Sigma,\Delta,Q_{0},F)$ where $Q$ is the set of states, $\Sigma$ is the alphabet, $\Delta:Q\times\Sigma\to 2^{Q}$ is the transition function, $Q_{0}\subseteq Q$ is the set of initial states, and $F\subseteq Q$ is the set of accepting states. Likewise, a Deterministic Finite Automaton (DFA) is a tuple $(Q,\Sigma,\delta,q_{0},F)$, defined analogously except that $\delta:Q\times\Sigma\to Q$ is the transition function and there is only a single initial state $q_{0}$.

We consider the NFA $M$ and the DFA $M^{\prime}$ described in Algorithm 1. The NFA $M$ has states indexed by strings in $\{A,C,G,T\}^{k-1}$ and alphabet $\Sigma=\{0,1,\ldots,b-1\}$, the current levels. Given a state $(s_{0}\cdots s_{k-2})$ and an input current level $i\in\Sigma$, the NFA $M$ can transition to any other state of the form $(s_{1}\cdots s_{k-1})$ so that $f(s_{0}s_{1}\cdots s_{k-1})=i$. All states are accepting states. By construction, a sequence of current readings $c\in\Sigma^{n-k+1}$ can be an output $f^{*}(s_{0}s_{1}\cdots s_{n-1})$ of the nanopore sequencer if and only if $c$ is accepted by $M$. The DFA $M^{\prime}$ accepts exactly the same strings as $M$, and is obtained using the classic subset construction [18], such that each state of the DFA corresponds to a subset of the states of the NFA.

The transfer matrix method is a method for obtaining a generating function $g_{f}(z)$ for the number of states accepted by a given DFA. In more detail, given the transfer matrix $T$ for the automaton (so that $T_{i,j}$ is the number of transitions from state $i$ to state $j$), the transfer matrix method gives an expression for a function $g_{f}(z)$ (in terms of the matrix $T$), so that

so that $N_{m}$ is the number of strings of length $m$ accepted by the finite automaton. We will use the notation $[z^{m}]g_{f}(z)$ to denote the coefficient $N_{m}$ on $z^{m}$.

In our case, the number of strings of length $m$ accepted by our finite automaton will be the number of current readouts of length $m=n-k+1$ that can be generated by some DNA strand of length $n$.

From $g_{f}(z)$, we can derive the asymptotic behavior of the number of possible current readouts required to determine $C_{f}$. As shown in the proof below, it is related to the smallest positive singularity of $g_{f}(z)$.

The above approach for computing $C_{f}$ exactly given a mapping $f$ is only practical for small window sizes $k$. However, we can derive some general bounds that apply to any mapping $f$. In particular, we are interested in the worst-case mapping $f$ (e.g., the one that attains $\min_{f}C_{f}$), as well as the best-case mapping $f$ (e.g., the one that attains $\max_{f}C_{f}$). We prove the following lemma.

Combining these bounds, we obtain the plot in Figure 3. One might wonder how tight the bounds shown in Figure 3 are, and also about where $C_{f}$ lies for a “typical” mapping function $f$. Using Algorithm 1, we are able to exactly calculate the maximum, minimum, and average $C_{f}$ for $k=2$, restricted to balanced mappings. These empirical results are shown in Figure 3.

The lower bound is on $\min_{f}C_{f}$ curious. Our theoretical bounds on $\min_{f}C_{f}$ are not tight, but the $k=2$ results suggest there may be some “bumpiness” in the true bound. Also based on the $k=2$ results, it appears that random (balanced) mappings are closer to the best-case scenario than the worst-case. Section 5.2 illustrates one coding scheme that takes advantage of a property shared by most random mappings for $b=2$, and we may hope there exist more general coding schemes that perform well for random mappings.

In the proof of Lemma 5, we observed that given a mapping $f$ with window size $k$ and $b$ distinct current levels, we can achieve a rate of $\frac{\log(b)}{k}$ by coding only on the $0^{th},k^{th},2k^{th},$ etc. current readings. Here, we propose two generalizations of this trivial scheme to improve the rate with additional preprocessing based on the mapping $f$.

We have initiated a theoretical study of coding for a highly abstracted version of the nanopore sequencer for DNA storage. We have provided algorithms and bounds for understanding the capacity, and we have given efficient coding schemes. However, we view our work as the tip of an iceberg. First, even for this abstracted model, much remains open. Can one derive better bounds on the capacity, or compute it efficiently for, say, $k=6$? Second, we hope that our insights will generalize to more practical models, including with substitution and synchronization errors.