Concatenated Codes for Recovery From Multiple Reads of DNA Sequences

Decoding sequences that stem from multiple transmissions of a codeword over an insertion, deletion, and substitution channel is a critical component of efficient deoxyribonucleic acid (DNA) data storage systems. In this paper, we consider a concatenated coding scheme with an outer low-density parity-check code and either an inner convolutional code or a block code. We propose two new decoding algorithms for inference from multiple received sequences, both combining the inner code and channel to a joint hidden Markov model to infer symbolwise a posteriori probabilities (APPs). The first decoder computes the exact APPs by jointly decoding the received sequences, whereas the second decoder approximates the APPs by combining the results of separately decoded received sequences. Using the proposed algorithms, we evaluate the performance of decoding multiple received sequences by means of achievable information rates and Monte-Carlo simulations. We show significant performance gains compared to a single received sequence.


I. INTRODUCTION
Coding-theoretic aspects of data storage in deoxyribonucleic acid (DNA) have recently gained significant attention. This new line of research has been triggered by a manifold of successful experiments [1]- [3] that demonstrated the viability of DNA as a reliable medium for archival data storage. Following the pioneering experiments, a variety of related information-theoretic problems with different foci have been identified. One of these problems is the reliable transmission over channels that introduce insertions, deletions, and substitutions (IDSs) [4], motivated by the fact that the sequencing process often introduces errors in the form of IDSs.
There is a variety of literature dedicated to the study of error-correcting codes for the IDS channel. While the list of related papers is long, the line of research most relevant to our work is the one inspired by Davey and MacKay [5]. In their paper, they propose a class of concatenated codes together with a decoding algorithm that represents the inner code and the channel by a hidden Markov model (HMM), which allows inferring a posteriori probabilities that are passed to the outer decoder. Following this work, the Markov process of convolutional codes was incorporated into that of the IDS channel, allowing to run standard decoding algorithms of convolutional codes for the IDS channel [6], [7]. The approach of [5] was later refined with an improved decoding algorithm [8]  using inner marker codes [9]. Recently, it was shown that synchronization is also possible without an inner code using a standalone low-density parity-check (LDPC) code [10].
The vast majority of previous studies has focused on the case of error correction after a single transmission over an IDS channel. In DNA-based storage, however, data is typically synthesized into many short DNA strands, and each of these strands is replicated thousands of times. Therefore, when accessing the stored data via DNA sequencing, multiple copies of each original strand are read. With this in mind, in this paper we deviate from previous work and consider error-correcting codes for multiple transmissions of a coded sequence over a channel that is impaired by IDS errors. In particular, we propose two new decoding algorithms for multiple received sequences and show that significant gains can be achieved as compared to a single transmission. The first decoding algorithm, mainly a benchmark for the second decoder, is an optimal symbolwise a posteriori probability (APP) decoder. While this algorithm has exponential complexity in the number of received sequences, we show that it is possible to achieve competitive performance with a second, practical decoder that has only linear complexity. For our results we employ a new concatenated coding scheme with an outer LDPC code and an inner convolutional code and compare with the watermark codes from [5]. To the best of our knowledge, this is the first work to propose decoding algorithms for the case of coded multiple transmission and to evaluate coding schemes as well as achievable information rates (AIRs) over this channel.
Related to our work is the study of sequence reconstruction, introduced by Levenshtein [11]. Here, a sequence is repeatedly transmitted over an erroneous adversarial channel and the goal is to analytically quantify, as a function of the sequence length, the number of sequences that are required to guarantee correct reconstruction of the transmitted sequence with zero error probability. Levenshtein's initial work on insertion and deletion errors has recently been extended by additionally allowing substitution errors [12].
Similar in spirit to our work, yet different in methods and objectives, is the research of trace reconstruction over deletion channels [13], [14], which is the probabilistic variant of the sequence reconstruction problem. While these works mainly derive asymptotic results, here we discuss a fixed number of short, finite-length received sequences over a channel that additionally allows insertions and substitutions. More recently, the trace reconstruction problem has also been formulated Delete: y ← y

Transmit
No error: for a fixed number of sequences with a larger focus on algorithmic aspects [15], [16]. However, these works consider only uncoded sequences, while we are interested in the case of coded transmission. Finally, in [17] Varshamov-Tenegolts codes have been proposed for error correction over multiple channels that introduce a fixed number of deletions.

A. Channel Model
We consider a state-based representation of the IDS process [5], [18] as depicted in Fig. 1. Let x = (x 1 , . . . , x N ), x i ∈ Σ q = {0, 1, . . . , q − 1} be the sequence to be transmitted over the channel. The transmitted symbols x i are successively queued for transmission and the received sequence y is assembled by the following procedure. When x i enters the queue, the channel enters state x i and there are three different possible events. First, a uniformly random symbol a ∈ Σ q is inserted to the received sequence y with probability p I . In this case, x i remains in the queue and the channel returns to state x i . Second, the symbol x i is deleted with probability p D without appending it to y and the next symbol x i+1 is enqueued. Third, x i is transmitted with probability p T = 1 − p I − p D . In this case, the symbol is substituted with a random symbol a * = x i with probability p S and the next transmit symbol x i+1 is enqueued. When the last transmit symbol x N leaves the queue, the procedure finishes and the channel outputs y. Note that y = (y 1 , . . . , y N ) has random length N , depending on the channel realization.

B. Forward Error-Correction Scheme
In this work, we will use a concatenation of an outer code with an inner code as our forward error-correction scheme, following [5]. The role of the inner code is to maintain synchronization and provide reliable soft information to the outer code, while the role of the outer code is to use this soft information to perform an accurate estimate of the transmitted information. In particular, given that the inner code provides reliable synchronization information, the outer code corrects the remaining errors on the synchronized sequence. In our scheme, the inner code is either a convolutional code as used in [19], or a watermark code as used in [5]. Note that the watermark code can be thought of as a convolutional code with memory equal to zero, and we will therefore restrict some passages to convolutional inner codes only. The outer code is either a binary or a nonbinary LDPC code.
The resulting codeword w = (w 1 , . . . , w No ), w i ∈ F 2 k is then encoded via either an [n, k, m] q inner convolutional code of blocklength n, binary input length k, memory m, and output alphabet Σ q or an [n, k] q inner block code. With this concatenation, each outer codeword symbol w i ∈ F 2 k is mapped to n symbols over Σ q . We terminate the convolutional code, resulting in the codeword The inner code consequently has rate R I = Nok N log q , and thus we transmit at an overall rate R = R o R I = Kk N log q . Before entering the IDS channel, the codeword v is offset by a pseudo-random sequence resulting in the final output x = (x 1 , . . . , x N ), x i ∈ Σ q . This offset sequence is known to the inner decoder, and is used to maintain synchronization at high probabilities of insertions or deletions. According to the nomenclature in [5] we call the combination of the inner block code and offset watermark code.
The output sequence x is then transmitted M times independently over an IDS channel resulting in received sequences y 1 , . . . , y M corresponding to M reads of the original strand. The inner decoder uses these received sequences to infer likelihoods for the symbols in w. These likelihoods are then fed to the outer decoder, which decides on the decoded sequenceû. Furthermore, we can also iterate between the inner and outer decoder, exchanging extrinsic information between them. The communication system is depicted in Fig. 2.

III. SYMBOLWISE A POSTERIORI DECODING OF INSERTION/DELETION/SUBSTITUTION CHANNELS
One of the main challenges for decoding from IDS channels is the random memory associated with the insertion and deletion process. This is visualized by the receiver's inability to directly identify the origin of a received symbol y i . Since insertions or deletions before symbol y i might have moved the symbol right or left in the received sequence, y i could be the result of transmitting a symbol x i with i = i. In the following, we describe how it is still possible to infer a posteriori likelihoods from this channel using a hidden Markov representation of the channel [5], [7], [8]. We present and discuss this HMM for M = 1 first and extend it to M > 1 multiple received sequences afterward. Here we restrict the discussion to convolutional inner codes, as the watermark code can be viewed as a convolutional code with m = 0.

A. Single Received Sequence
The APPs of the outer codeword symbols are given by The standard approach to compute p(y, w i ) for convolutionally encoded transmit sequences is to de-marginalize the distribution p(y, w i ) with respect to the memory states of the convolutional code that correspond to the information symbol w i and then exploit the Markov property of the convolutional code. That is, given the memory state at time i, the output after time i becomes independent of earlier states and the preceding output symbols. However, for the case of the IDS channel, this Markov property is violated through the memory of the channel. Therefore, it is useful to add a new hidden state variable, the so-called drift [5]. The key property of the drift is that, by its inclusion as a new state variable inside the Markov process, the Markov property is restored. This is because the drift sequence itself forms a Markov chain and, conditioned on the memory state and drift at time i, the output of the channel after time i becomes independent of previous memory states and drifts. Denoting the state variables of the convolutional code by s i , i = 0, 1, . . . , (N o +m), and introducing the joint state variable σ i = (s i , d i ), we obtain by slight abuse of notation p(y, w i ) = (σ,σ ):wi p(y, σ, σ ) , where the summation is over all state transition pairs (σ, σ ) that correspond to the information symbol w i . Further, due to the Markov property, we can expand the joint probability to p(y,σ,σ ) = p y (i−1)n+d 1 , σ p y in+d (i−1)n+d+1 , σ σ p y N in+d +1 σ . Abbreviating the above terms by α i−1 (σ), γ i (σ, σ ), and β i (σ ) in order of appearance, one can deduce the forward and backward recursions with appropriate termination conditions incorporating the termination of the convolutional code and drift states. The branch metric γ i (σ, σ ) = p(w i )p(y in+d (i−1)n+d+1 , d d, s, s ) can be efficiently computed using a lattice structure, see, e.g., [7], [20] 1 using the transmit sequence corresponding to the state transition s → s . Having these expressions in hand it is possible to compute the APPs p(w i |y) using the BCJR algorithm [21]. Note that while the introduction of the drift state variable allows to exactly compute the APPs, it comes at the cost of additional complexity in terms of a larger number of states and state transitions, as discussed in Section III-C.

B. Multiple Received Sequences
In this work we consider the following two approaches of computing the a posteriori distribution p(w i |y 1 , . . . , y M ) for multiple sequences.
Joint decoding: Similar to the approach for a single received sequence, it is possible to introduce a drift state sequence The forward, backward recursion, and the terminal conditions are according to the single sequence case.
Separate decoding: Instead of computing the a posteriori likelihoods precisely, we now propose a different approach that allows to approximate the joint a posteriori likelihood using the marginal likelihoods p(w i |y j ) as The approximation is motivated by the case of memoryless channels, where the proportionality holds exactly. The individual APPs p(w i |y j ) can be computed as in the case of a single received sequence by means of a BCJR decoder. While the approximation comes at reduced complexity as compared to the joint decoding approach, see Section III-C, it also results in a loss of performance, which we will quantify in Section IV. Note that it is also possible to separate the joint a posteriori likelihoods at other stages, such as directly between the channel and inner code or after decoding the LDPC code. However, we observed that for the considered parameters the proposed stage yields the best results.

C. Complexity Considerations
The number of operations performed by the inner BCJR decoder is proportional to the total number of edges of the trellis that is inferred from the HMM. In order to reduce its computational complexity, the drift states are usually limited to a fixed interval d i ∈ [d min , d max ] and the maximum number of insertions per symbol is limited to I max [5]. Denote by ∆ = d max − d min + 1 the total number of drift states and by δ = n(I max + 1) + 1 the number of possible drift transitions. The complexity, measured by the total number of edges in the trellis, of joint decoding is then given by The [1, 1, 2] 4 convolutional code is a binary convolutional code with generator polynomials (1 + D 2 , 1 + D + D 2 ), where each two bits are combined to one DNA symbol. The [3,3] 4 watermark code is constructed as in [5].
while the complexity of separate decoding is simply given by M times the complexity of decoding a single sequence, Note that here we neglect the effect of trellis termination on the number of edges. Typical values of ∆δ for the DNA storage application can easily be as high as ∆δ ≈ 1000. This means that decoding the inner code is the main contributor to the decoding complexity and we thus neglect the complexity of the outer decoder in the following. It further becomes evident that the joint decoding approach only works for a very small number of sequences M due to its exponential complexity.
IV. SIMULATION RESULTS In our simulations, we compare the two following set-ups. First, an inner convolutional code with an outer WiMAX binary LDPC code (CC + B-LDPC) and, second, an inner watermark code with an outer WiMAX-like nonbinary LDPC code (WM + NB-LDPC). Their performance is compared over a single and multiple IDS channels. Recall that in both schemes a pseudo-random offset sequence is used to improve the synchronization capability.

A. Simulation Parameters
Our simulations are done over an alphabet of size q = 4, which equates to the four bases {A, C, G, T} used in DNA and we set p S = 0 and p I = p D throughout. The code parameters are summarized in Table I. The outer LDPC code is decoded via belief propagation with a maximum number of 100 iterations. Regarding the inner code, we set I max = 2 and the limit of the drift random variable in decoding is set to be five times the standard deviation of the final drift at position N , i.e., d max = −d min = 5 √ N p D 1−p D . For multiple transmissions, we use M = 2, 5, 10, and 20. Fig. 3 shows the frame error rate (FER) of the CC + B-LDPC scheme for single and multiple transmission over an IDS channel, while Fig. 4 shows that of the WM + NB-LDPC scheme. Interestingly, when transmitting over a single IDS channel, the CC + B-LDPC scheme outperforms the WM + NB-LDPC scheme. However, for transmission over multiple IDS channels, this behavior reverses, and the WM + NB-LDPC scheme has superior performance for a larger number of sequences. While this result seems surprising at the first glance, there is a natural explanation for this phenomenon as follows. First notice that with the inner convolutional code, bitlevel APPs p(w i |y j ) with w i ∈ F 2 are computed, while the watermark code works over nonbinary inputs w i ∈ F 2 3 . Now, for larger error probabilities p I , p D , the memory of the channel  becomes more apparent and thus also the a posteriori distribution has increased correlation between different input symbols w i . While this correlation is marginalized for the binary APPs, it is better captured by the nonbinary input symbols and thus the WM + NB-LDPC scheme achieves better performance for larger number of sequences or higher error probabilities. In terms of complexity, both schemes have the same ∆ and M , but different n, k, m, and δ. Comparing both complexities for equal blocklengths, we get C conv.

B. Discussion of Numerical Results
sep /C watermark sep ≈ 1.2. Thus, with these parameters the convolutional code is slightly more complex than the watermark code. Since the complexity increases exponentially with the memory m, larger memory convolutional codes are expected to become even more complex as compared to the watermark code. The same conclusion can be made when multiple channel transmissions are combined with inner-outer iterative decoding, since doing M iterations has the same complexity as separate decoding of M received sequences. Fig. 5 shows the FER of both schemes when iterating times between the inner and outer decoder. We fix the total complexity such that the product of the number of iterations and the number of sequences is equal to 20. It is observed that the iterative scheme benefits more from a larger number of received sequences, however in practice this comes at the cost of additional sequencing costs. Finally, AIRs with an inner convolutional code are displayed in Fig. 6. In agreement with our previous results, the achievable rates are significantly larger for multiple transmissions as compared to a single transmission. Not surprisingly, the joint decoder achieves higher rates. It exhibits a similar performance with M = 2 transmissions as compared to the separate decoder with M = 5 sequences.

V. FUTURE WORK
In view of the huge number of trellis states that the joint decoding approach entails, an immediate idea is to consider suboptimal decoding of the inner code, such as sequential  decoding. There are however some obstacles that have to be overcome, including a huge number of edges per node and an appropriate metric. This will be addressed in future work.