Keywords: DNA, information storage, error-correcting code, indel, Reed–Solomon

HEDGES Theoretical Design.

Fig. 1 shows the data flow for HEDGES encoding and decoding algorithms. In the terminology of coding theory, HEDGES is an infinite-constraint-length convolutional code (a “tree code”) incorporating some features specific to the DNA channel. It is decoded via a stack algorithm that assigns costs to both indels and substitutions. The decoding algorithm succeeds probabilistically, with the ability to signal success or failure. Decoding failures are then regarded as erasures (unknown bits or bytes) and can be corrected in the outer code [i.e., an RS(255,223) code]. Alternatively, the error strand can be discarded and resequenced.

The simplest case is a half-rate code (1 bit encoded per nucleotide) with no constraints on the output DNA sequence, shown inFig. 1D (seeSI Appendix, Fig. S1 for full diagram). The basic plan is a variant of a centuries-old “text auto-key encoding” cryptographic technique (16) (see also Wikipedia, “Autokey cipher”). We generate a stream of pseudorandom characters$K_i\in \left\{\mathrm{0,1,2,3}\right\}\equiv \left\{A,C,G,T\right\}$, where each$K_i$ depends deterministically, via a hash function, on a fixed number of previous message bits$\left\{b_j\right\}$, on the current bit position index$i$ of the current message bit$b_i$, and on the strand ID. We then emit a character$C_i$ in$\left\{A,C,G,T\right\}$ with$C_i=K_i+b_i$, the addition performed modulo 4. Redundancy for error correction occurs because, at each$i$,$C_i$ can take on only two out of four values, because$b_i\in \left\{\mathrm{0,1}\right\}$. The output DNA is always completely (pseudo)random because the hash is pseudorandom. Thus, the DNA is, in this sense, independent of the encoded message. We note that nothing limits this scheme to the four “natural” nucleotides; generalization to an increased number is straightforward.

The decoding algorithm sequentially guesses message bits (Fig. 1F). Each guessed bit$b_i$—along with its guessed position$i$—allows the algorithm to predict (via the forward encoding algorithm) a nucleotide value$C_i$. If$C_i$ agrees with the observed value, a reward is assigned to that guess. If$C_i$ disagrees, the guess is penalized as appropriate for a substitution. If it disagrees, but would agree at a different assumed position, plus or minus one, that distinct guess is assigned a penalty appropriate for an insertion or deletion. A heap structure keeps systematic track of all currently viable chains of guesses. To prevent the heap from growing exponentially, only chains of guesses that have close to the best total score are extended preferentially, using a variant of the$A^{*}$ algorithm (17).

In summary, the algorithm encodes information as a stream of nucleotides such that any single decoding error in either nucleotide identity or nucleotide position will “poison” the downstream predictions. Thus, on decoding, there will be only one good-scoring chain of guesses—the correct one. In the unlikely case that the heap grows larger than a preset size$H_{\text{limit}}$ (e.g.,$10^6$), we declare a decode failure and mark the remaining part of the strand as an erasure. Similarly, by including the strand ID in the hash function input, any strand with incorrectly decoded strand ID will be “poisoned” for the full length of the message and fail to decode. This results in effectively zero strand ordering errors, the most expensive type of simple error in interleaved encoding schemes.

Testing In Silico.

For in silico testing, we assumed equal rates for substitutions, insertions, and deletions with total error probability per nucleotide$P_{\text{err}}=P_{\text{sub}}+P_{\text{ins}}+P_{\text{del}}$. Any other distribution can then be conservatively bounded by the choice$P_{\text{err}}=3max\left(P_{\text{sub}},P_{\text{ins}},P_{\text{del}}\right)$. For simplicity, we assumed that errors occurred at random positions in the DNA strand. The experimental tests in vitro had no such assumption (see below).

The HEDGES algorithm was implemented in C++ for speed, with a Python-callable interface for encoding/decoding single strands. As an initial validation of programming accuracy and interface design, we used an outer-code concatenated design with packets of 255 strands of length 300 protected by RS(255,223). A detailed description of the encoding design is provided inSI Appendix,Supplementary Text; all corresponding computer programs are provided inSI Appendix and available via GitHub. For every code rate$r$ in$\left\{0.166,0.250,0.333,0.500,0.600,0.750\right\}$, and each assumed$P_{\text{err}}$ in$\left\{0.01,0.02,0.03,0.05,0.07,0.10,0.15\right\}$, we encoded 10 packets ($\sim 10^6$ nucleotides), pooled all of the strands, and duplicated them as if sequencing to depth 5, meaning that each strand was duplicated a Poisson random number of times, with mean 5. We verified that the pooled and duplicated strands could be decoded, the packets recovered and ordered, RS applied, and the messages recovered without error—but only up to some maximum tolerable error rate$P_{\text{err}}$ that increased with decreasing code rate. That is, we “tested to failure” on$P_{\text{err}}$. For this test, the maximum heap size was set to an intentionally small value,$H_{\text{limit}}=10^5$, to increase the number of decode failures and thus stress the program. The results of this test were as expected and gave us confidence to proceed with in vitro and larger-scale in silico testing (below).

We next needed to construct a statistical error model that could be extrapolated to the petabyte or exabyte scale. For this model, we needed to know the rate of bit errors and byte errors in HEDGES output (for each code rate$r$ as a function of$P_{\text{err}}$). Because decode errors tended to occur in bursts, the rate of byte errors was less than the 8$\times$ expected for independent bit errors. We also needed to know the probability per strand of a decode failure, and the distribution of such failures along the strands. For each pair$\left(r,P_{\text{err}}\right)$, we simulated 12,000 strands of length 10,000 (about$10^8$ nucleotides), now with$H_{\text{limit}}=10^6$. Decode failures occurred approximately uniformly along the strands, consistent with the hypothesis that the heap decode “loses its way” only on rare, local, random sequences (SI Appendix, Fig. S2). Decode failures are thus characterized by a single value, the mean run length to failure in a Poisson model.Fig. 2A shows the byte error rate as a function of code rate, whileSI Appendix, Fig. S2 gives full details on observed bit and byte error rates, and mean runs to decode failures. Byte error rates were typically 3 to 5 (not 8) times the bit error rates.

Using these byte error rates, we then modeled HEDGES in an overall concatenated ECC design.Fig. 2B shows the average number of message bytes that could be decoded before encountering an uncorrectable error using the concatenated design previously described (seeMethods for details). A broad set of code rates (green) are suitable for gigabyte- to exabyte-scale DNA storage. HEDGES decode failures in this region occur every$10^4$ to$10^5$ nucleotides but get corrected by the RS outer code or, optionally, by using another exemplar of the strand (seeDiscussion). Similar simulations with the imposed output constraints of no homopolymer runs (e.g., GGGG or AAAA) greater than four, and$4\le \text{CG}\le 8$ in any window of 12 nucleotides, are shown inSI Appendix, Fig. S3, and are not substantially different fromFig. 2. We also modeled the effects of sequencing constraints more generally (SI Appendix,Supplementary Text and Fig. S4), and the combined model and simulation results indicate that the most common sequencing constraints have minimal impact on HEDGES. In sum, in silico simulations indicate that HEDGES is capable of error-free decoding of exabyte-scale messages.

Testing In Vitro.

We next tested real-world ECC performance on a pooled sample of 5,865 synthetic 300-base pair DNA strands that were exposed to accelerated aging or enzymatic mutagenesis. Of these, 18 packets of 255 strands were HEDGES inner-encoded (with subsets at each of six code rates) and then RS(255,223) outer-encoded across strands. Five packets, totaling 1,275 strands, were encoded with an unrelated error correction algorithm (18), but also served as a negative control on identifying and sequencing HEDGES strands into packets. Each HEDGES strand consisted of$3^{\prime }$ and$5^{\prime }$ primers of length 23 nucleotides (seeMethods) flanking a 254-nucleotide DNA payload. When decoded into bytes, each payload comprised a 1-byte packet number, a 1-byte sequence number (these “salt-protected” on encryption; seeMethods), a message payload whose length depended on the code rate, and a 2-byte runout. The sample was PCR amplified and prepared for Illumina-based sequencing. Additionally, we degraded the DNA separately via error-prone PCR mutagenesis or by incubation at high temperature (seeMethods). Sequencing was done to a mean depth of$\sim$50.

We performed two kinds of tests of the decoding algorithm, with and without knowledge of the encoded message. “Type A” tests assumed knowledge of the 5,865 strand sequences and could be used to characterize the nature of end-to-end DNA error rates. “Type B” tests were blind decodings of the sequenced data, with knowledge only that the pooled DNA contained HEDGES-encoded data in the specified format.

In our Type A tests, 10 to 15% of sequenced strands could not be uniquely identified with any known input strand, even using quite robust N-gram methods, and even for unmutagenized aliquots. This may be the result of the low concentration, or of contamination at some stage; but it also added to the challenge for the blind type B tests.

For strands whose progenitor sequence could be identified,Table 1 shows measured rates of substitution, insertion, and deletion errors. Notably, only the highest mutagenesis kit protocol produced a substantial increase in DNA errors. Data in ref.3 estimate DNA degradation over a wide range of timescales and temperatures, suggesting that 50 °C incubation for 8 h should have produced significant mutagenesis. We did not find this, however. So, for further analysis here, we consider only the untreated and high-mutagenesis datasets.

Table 2 shows the results for decoding strands that were identified as belonging to packets of each code rate. Approximately 3% of the strands failed to decode even at small code rates where, in simulation, there were many fewer such failures. Identification of these strands was ambiguous and may stem from PCR mispriming, oligonucleotide misdimerization, and other next generation sequencing (NGS) library preparation artifacts that can vary from batch to batch. Indeed, for lower code rates (where the ECC was relatively unstressed), strand decode failure rates were slightly higher for the untreated case than for the high-mutagenesis case, presumably due to batch-to-batch variation in the number of such artifacts.

For this reason, the values for the mean run to an uncorrectable error inTable 2 are calculated assuming a strategy of rejecting failed decodes, rather than counting them as erasures. We adopted just this rejection strategy in our blind (type B) decoding; the input data were several$\times 10^5$ total reads of the synthesized 5,865 strands (of which 4,590 were message bearing) plus contamination. We shuffled the reads into random order, then attempted HEDGES decoding one read at a time, populating the 18 expected packets of 255 strands with successful decodes and attempting the outer RS error correction when the number of erasures (missing strands) was small enough (seeMethods for details).

As expected based on the results ofTable 2, we achieved error-free decodes of all packets, except in the case of two packets with high mutagenesis at the highest code rate 0.750. With no mutagenesis, 24,000 total reads were required for all 18 packets. With high mutagenesis, 22,000 reads were required for 16 packets, while the undecodable two continued to fail indefinitely. In the successful cases, the number of reads corresponded to about depth 3 on message-bearing message strands. This depth was required merely to sufficiently populate the packets for the outer code to operate due to random strand sampling, not because of any property of HEDGES as the inner code.

HEDGES Encoding in the Half-Rate Case.

Given a message stream of bits

(“the message” or “bits”), we want to emit a stream of DNA characters

(“the codestream” or “characters”). We first review the case of a half-rate code, where we emit exactly one$C_i$ (2 bits of output) for each$b_i$ (1 bit of input). Then we will generalize to codes at other rates$r$ (message bits per codestream bit),$0<r<1$, so that the streams$b_i$ and$C_i$ are not then in lockstep, and$M\hspace{0.17em}\ne \hspace{0.17em}N$. One should think of$N$ as being on the order of$10^2$ to$10^4$, the maximum length of a single DNA strand that can be cheaply synthesized today or in the foreseeable future. We want to be able to decode, without residual errors, a received codestream$C^{\prime }$ that differs from$C$ by substitutions (errors), insertions, and deletions (collectively “indels”). Indels are silent: Their positions in the codestream$C^{\prime }$ are not known to the receiver.

We generate a keystream of characters$K_i\in \left\{\mathrm{0,1,2,3}\right\}$, where each$K_i$ depends pseudorandomly (but deterministically by a hash function) on some number of previous message bits$b_j$ (with$j<i$), and also directly on the bit position index$i$ and the strand ID,

where$S_i$ is$s$ bits of salt (strand ID),$I_i$ is the lower$q$ bits of the bit index,$B_i$ is the previous$p$ bits, and$F$ is a chosen hash function. (SeeSI Appendix,Supplementary Text for initialization.) We then emit a codestream character

the addition performed modulo 4.

The redundancy necessary for error correction comes from the fact that$b_i$ takes on only two values, while$K_i$ and$C_i$ can have four values. This generates (only) 1 bit of redundancy per character, that is, can be acausally valid by chance half of the time. However, the dependence of$K_i$ on many previous message bits ties any given message bit to many future bits of redundancy. Similarly, the dependence of$K_i$ on$i$ ties every bit to its position index, so that insertions can be identified and removed, and deleted values can be restored.

Some further details about the encoding algorithm are given inSI Appendix,Supplementary Text.

HEDGES Decoding Algorithm.

For simplicity, assume that error rates are “small,” so that “most” DNA bases are received as they were intended. (We saw, inResults, that DNA character error rates up to$\sim 5$ to 10% are tolerable.) Suppose we have correctly decoded and synchronized the message through bit$b_{i-1}$ and now want to know bit$b_i$. Guessing the two possibilities,$\left\{\mathrm{0,1}\right\}$, we useEq.3 to predict two possibilities for the character$C_i$. In the absence of an error, only one of these is guaranteed to agree with the observed character$C_i^{\prime }$. We assign, to a guess that generates disagreement with$C_i^{\prime }$, a penalty score equal (conceptually) to the negative log probability of observing a substitution error. In other words, a wrong guess might actually be right, but only if a substitution has occurred. If neither guess produces the correct$C_i$, then both are assigned the substitution penalty.

We have not yet accounted for the possibility of insertions and deletions, however. In fact, there are more than the above two possible guesses. We must guess not just$b_i\in \left\{\mathrm{0,1}\right\}$, but also a “skew”$\mathrm{\Delta }\in \left\{\dots ,-\mathrm{1,0,1},\dots \right\}$ that tells us whether, in comparing$C$ to$C^{\prime }$, we should skip characters ($\mathrm{\Delta }>0$) because of insertions, or posit missing characters ($\mathrm{\Delta }<0$) because of deletions (in which case, there is no comparison to be done). As a practical simplification, we consider only$\mathrm{\Delta }\in \left\{-\mathrm{1,0,1}\right\}$, requiring multiple indels to resolve as concatenated single indels. Then there are six guesses for$\left(b_i,\mathrm{\Delta }\right)\in \left\{\mathrm{0,1}\right\}\otimes \left\{-\mathrm{1,0,1}\right\}$. Each can be scored by an appropriate log probability penalty for any implied substitution, insertion, or deletion.

Log probability penalties accumulate additively along any chain of guesses. In the causal case of a chain of all-correct guesses, we accumulate penalties only in the (relatively rare) case of actual errors. However, because of the way that the key$K_i$ (Eq.3) is constructed, a single wrong guess for either$b_i$,$i$, or$\mathrm{\Delta }$ throws us into the acausal case where 3/4 of subsequent comparisons of computed$C$ (at some bit position index$i$) to observed$C^{\prime }$ (at some index$k$) will not agree—thus penalties will accumulate rapidly. The decoding problem, conceptually a maximum likelihood search, thus reduces to a shortest-path search in a tree with branching factor 6, but with the saving grace that the correct path will be much shorter than any deviation from it.

The rate of decode errors rises in the last several bytes of message, because some incorrect chains don’t have time to accumulate bad scores. To counter this, we pad each strand with (typically) two “runout bytes” of message zeros at encode, and ignore them at decode. The need for runout bytes makes the HEDGES algorithm inefficient (and thus unsuitable) for an application needing very short DNA strands (e.g., tens rather than hundreds of nucleotides).

Further details about the decoding algorithm are given inSI Appendix,Supplementary Text.

Use of Salt to Protect Critical Message.

InEq.3 andFig. 1D, we allowed for some number of bits of known salt$S_i$ when message bit$b_i$ is encoded. The use of salt is optional, but is desirable in an overall interleaved design where sequenced strands from a pool need to be correctly ordered at decode time (SI Appendix,Supplementary Text). This is generally the case when the outer code is interleaved across strands. To the outer code decoder, each incorrectly ordered strand is equivalent to a full strand length of random errors, so it is very important to protect strand ID message bits that determine the strand ordering for outer decoding. Here is how salt is enabling of extra protection: Suppose we want to protect an initial$s$ message bits. Then define, recursively, the salt by

Most errors in the first$s$ bits will be corrected as usual by the shortest-path heap search. But any residual error that gets through will “poison” the salt for the entire rest of the strand, rendering it undecodable. In effect, we convert an error in the protected bits into an erasure of the whole strand. This may seem drastic, but it is just what we want: A strand with incorrect serial number (and hence incorrect ordering among other strands) would look like a strand of errors to the outer code; an erased strand is equivalent to only half as many errors.

Code Rates Other than One-Half.

A simple modification of the encode and decode algorithms described above allows for code rates other than one-half. Take the input bitstream of expression1 and partition it into a stream of values$v_k$ with variable numbers of bits in the range 0 to 2, according to a repetitive pattern like the ones shown inTable 3. See alsoSI Appendix, Fig. S1.

Here are two examples showing how to interpret the entries inTable 3 (with adjacency denoting 2-bit values in${\mathbb{Z}}_4$):

Eq.3 for encoding now becomes

where$V_i$ is composed of concatenated previous variable bits. Pattern values of 0 provide 1 bit of additional redundancy check relative to the base case of code rate one-half, while pattern values of 2, encoding 2 bits per DNA character, provide 1 bit less (i.e., zero). By construction, the code rate is one-half the average of the integers in the pattern. The column in the table labeled$P_{\text{ok}}$ is a “greediness parameter” that mitigates the tendency of the heap to expand exponentially; seeSI Appendix,Supplementary Text for details of this.

Decoding follows exactly the same pattern. Guessing a 2-bit$v_i$ spawns 12 child hypotheses, while guessing a zero-bit$v_i$ spawns only 3.

HEDGES Parameters.

For encoding, the parameter choices are 1) the choice of code rate and variable bit pattern (as inTable 3), the default case being code rate 0.5; 2) the number$q>0$ of low-order bits of position index in the hash; 3) the number$p>0$ of previous message bits in the hash; 4) the number$s\ge 0$ of salt bits; and 5) the number$n\ge 0$ of initial message bits to be protected by salt. Aspects of the choices of, and trade-offs among, these parameters are further discussed inSI Appendix,Supplementary Text.

It is an important point that choosing the decode runtime parameters, for example,$H_{\text{limit}}$ or$P_{\text{ok}}$, is not an irrevocable choice. Given a DNA message, one can make multiple tries, varying the decode parameters adaptively until acceptable performance is achieved. Simply increasing$H_{\text{limit}}$ and retrying will often rescue failed decodes (SI Appendix,Supplementary Text and Fig. S6). One can evaluate success by running time and by the count of errors needing correction by the outer RS code. The parameter values that we suggest may be viewed as starting points.

Imposing DNA Output Sequence Constraints.

DNA synthesis and sequencing platforms have sequence-dependent error profiles. Imbalanced GC content and homopolymer runs are well known to be problematic, for example, leading to indel and substitution errors or even whole strand dropout errors in popular sequencers such as those from Illumina and Oxford Nanopore (19,20). Thus, many proposed ECCs impose constraints on GC content, homopolymer runs, or both. These typically involve one-off coding designs for each constraint, often reducing significantly the effective code rate (7). An important property of the HEDGES algorithm is that it can accommodate a large class of sequence constraints without additional, one-off, code design. Moreover, constraints may be imposed without decreasing the code rate (a seeming paradox that we explain below).

Consider the class of constraints that can be applied to a nucleotide sequence with only “past” information. That is, from the emitted sequence$\left[\dots ,C_{i-3},C_{i-2},C_{i-1}\right]$, we can determine whether the choice of$C_i$ is constrained and, if so, what is the set of its allowed values, here denoted$\left\{C_i^{*}\right\}$, whose size we denote$\#C_i^{*}$. We assume$0<\#C_i^{*}\le 4$, because a zero value would imply a constraint so severe that the strand cannot be continued at all. (One could adopt the convention of relaxing such constraints if they are sufficiently rare.)

The only required change in the HEDGES algorithm is to replaceEq.5 by

and output the thus-indexed character in the acceptable set$\left\{C_i^{*}\right\}$ (SI Appendix, Fig. S1).

Note thatEq.6 has the same code rate as [5]. How is this possible when constraints always act to reduce the number of possible output strings, and hence reduce the channel capacity? The answer is that [6] absorbs the reduced capacity completely into the error correction, not into the message. The extreme case is when$\#C_i^{*}=1$. Then, the emitted character does not depend on the message bits$v_i$, which become de facto erasures. At decode time, the missing bits are restored by virtue of the fact that they enter$V_i$ and thus the hash function for later emitted characters.

This is a versatile scheme. When there is a possibility of emitting an unacceptable homopolymer,$\#C_i^{*}$ will decrease from 4 to 3. When the GC count in a specified window is too large (or small),$\#C_i^{*}$ will decrease from 4 to 2. More-complex constraints are easily added. For example, Illumina sequencers have been reported to have high error rates following a GGC motif (4). To disallow GGC motifs, one might have imagined building a bespoke code around triplets of characters encoded as members of a large Galois field (3) and removing GGC, NGG, and GCN as possible output triples. Within HEDGES, GGC can be disallowed as an afterthought, simply by adding it to the list of forbidden outputs.

Statistical Model of System Designs with Outer Codes.

We take as input, fromSI Appendix, Fig. S5, the byte error rate$P_{\text{byte}}$ and mean run to decode failure$L_{\text{fail}}$. The outer code is characterized by its length$N_T=N_m+N_c$, where$N_m$ is the number of (payload) message bytes, and$N_c$ is the number of correction bytes. For RS(255,223), we have$N_T=255$ and$N_m=223$. The number of correctable byte errors is well known to be$N_c/2$, with erasures (bytes known to be missing) counting as half of an error.

We assume that the RS error correction is interleaved, that is, runs across strands. We further assume that, as in our example concatenated design, each RS input byte string samples uniformly along the length of the DNA strands and thus sees uncorrelated errors and erasures with their respective mean rates. We implemented this by applying the RS outer code “diagonally” across strands, so that strand ends (for example) are distributed over multiple RS decodes (Fig. 1C). The number of errors in each RS correction is thus a random Poisson variable${\mathcal{N}}_{\text{err}}$ with mean$N_T{P}_{\text{byte}}$.

Erasures can occur for two reasons: 1) If the DNA pool is sequenced to a mean depth$D$, then a fraction$p_{f1}=\exp \left(-D\right)$ strands (and therefore bytes) will occur zero times and be counted as erasures. 2) Along each strand of length$S$, the probability of a (first) decode failure at position$x$ is$\left(1/L_{\text{fail}}\right)\exp \left(-x/L_{\text{fail}}\right)$, in which case$S-x$ DNA positions are left undecoded. The fraction (in either DNA or bytes) of such erasures is thus

The number of erasures is thus a random Poisson variable${\mathcal{N}}_{\text{era}}$ with mean$N_T\left(p_{f1}+p_{f2}\right)$. The mean run (in message bytes) to an uncorrectable error is thus the reciprocal of a tail probability,

While the sum of two Poisson variables is Poisson, the factor one-half in front of${\mathcal{N}}_{\text{era}}$ spoils this, soEq.8 involves no simple distribution. Replacing the one-half by either zero or 1 does give Poisson distributions, however, that bound the desired result and turn out to be not too different numerically (in their power-of-ten exponents, which is all we care about). As an approximation, we thus may interpolate by setting the factor back to one-half and pretend that it is still Poisson.

Poisson tail probabilities can be written exactly as incomplete gamma functions (21). In our regime of interest, the upper and lower bounds (22),

nearly coincide, where$p_n=\text{Prob}\left(x=n\right)=\exp \left(-\lambda \right){\lambda }^n/n!$. We can use the upper bound to get a tight lower bound on$L_{\text{u.c.}}$.

The values inFig. 2 (red/yellow/green cells) result from performing the above calculations with$N_T=255$,$N_m=223$,$N_c=32$,$S=300$, and$D=10$.

Blind Decoding Tests.

We eliminated reads in which the flanking$5^{\prime }$ or$3^{\prime }$ primers were not found, or were separated by significantly more or less than the expected payload length of 254. Because any blind test must be ignorant of the message, we did not eliminate the 10 to 15% of strands that (we secretly knew) contained “junk” payloads that corresponded to no intended message.

Taking the reads one at a time, we first attempted decodes at all six designed code rates, but limiting the hypothesis budget to$H_{\text{limit}}=\mathrm{5,000}$. Whichever code rate got farthest in the message was declared as the putative actual code rate for that strand. At that code rate, we next attempted to decode 2 bytes (the packet number and strand sequence IDs) plus three runout bytes. If this decode succeeded, we checked by IDs to see whether this strand was already seen and decoded. If not, we performed a full decode with a hypothesis budget of$H_{\text{limit}}=10^6$. This strategy was designed to avoid unnecessary full decodes.

When the full decode succeeded, we accepted it as authoritative for that strand by packet and sequence ID, and performed no further full decodes for that strand when encountered. When it failed, we discarded the strand. That is, we did not attempt to rescue bytes up to the failure point.

After every 1,000 strands, we looked for packets that were sufficiently populated to possibly allow RS(255,223) outer error correction, counting missing strands as erasures. When the outer correction eventually succeeded, we marked that packet as done. This is a conservative strategy for testing, using only the first observed read of each strand and a near-maximal number of missing strands in each packet for the outer decode. In practice, after sequencing, one would use all available reads and potentially reduce the error rates even further.

DNA Library Prep and Sequencing.

A DNA library of 5,865 300-nt oligonucleotides was synthesized by Twist Biosciences. To introduce mutations, the DNA storage library was mutagenized using the Diversify PCR Random Mutagenesis Kit (Takara 630703). Samples were added with 480$\mathrm{\mu }$M${\mathrm{M}\mathrm{n}\mathrm{C}\mathrm{l}}_2$ and 40$\mathrm{\mu }$M 2′-deoxyguanosine 5′-triphosphate (dGTP), 640$\mathrm{\mu }$M${\mathrm{M}\mathrm{n}\mathrm{C}\mathrm{l}}_2$ and 200$\mathrm{\mu }$M dGTP, or 960$\mathrm{\mu }$M${\mathrm{M}\mathrm{n}\mathrm{C}\mathrm{l}}_2$ and 240$\mathrm{\mu }$M dGTP to achieve low, medium, and high levels of mutagenesis, respectively. Samples were then PCR amplified for 12 cycles using primers IF538 and IF539.

To mimic time-dependent aging of DNA, DNA storage library samples were incubated at 50 °C for a period of 2 h or 8 h in water (3). Samples were then PCR amplified for 12 cycles using primers IF538 and IF539.

PCR-amplified samples were purified using a PCR cleanup kit (NEB T1030).

Illumina sequencing libraries were prepared using the NEBNext Multiplex Oligos for Illumina (E7335) primer set. Libraries were sequenced on a paired-end MiSeq chip with 300 base pairs of read length. Final library preparation and sequencing was performed by the University of Texas Genomic Sequencing and Analysis Facility (GSAF).

• 1.Church G. M., Gao Y., Kosuri S., Next-generation digital information storage in DNA. Science337, 1628 (2012) [DOI] [PubMed] [Google Scholar]
• 2.Goldman N., et al. , Towards practical, high-capacity, low-maintenance information storage in synthesized DNA. Nature494, 77–80 (2013). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 3.Grass R. N., Heckel R., Pudda M., Paunescu D., Stark W. J., Robust chemical preservation of digital information on DNA in silica with error-correcting codes. Angew. Chem. Int. Ed.54, 2552–2555 (2015). [DOI] [PubMed] [Google Scholar]
• 4.Bornholt J., et al. , A DNA-based archival storage system. Comput. Architect. News44, 637–649 (2016). [Google Scholar]
• 5.Erlich Y., Zielinski D., DNA Fountain enables a robust and efficient storage architecture. Science255, 950–954 (2017). [DOI] [PubMed] [Google Scholar]
• 6.Yazdi S. M. H. T., Gabrys R., Milenkovic O., Portable and error-free DNA-based data storage. Nat. Sci. Rep.7, 5011 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 7.Organick L., Ang S. D., Chen Y.-J., Random access in large-scale DNA data storage. Nat. Biotechnol.36, 242–248 (2018). [DOI] [PubMed] [Google Scholar]
• 8.Ceze L., Nivala J., Strauss K., Molecular digital data storage using DNA. Nat. Rev. Genet.20, 456–466 (2019). [DOI] [PubMed] [Google Scholar]
• 9.MacWilliams F. J., Sloane N. J. A., The Theory of Error-Correcting Codes (North Holland, 1983). [Google Scholar]
• 10.Roth R. M., Introduction to Coding Theory (Cambridge University Press, 2006). [Google Scholar]
• 11.Moon T. K., Error Correction Coding: Mathematical Methods and Algorithms (Wiley, 2005). [Google Scholar]
• 12.Lin S., Costello D. J., Error Control Coding (Pearson, ed. 2, 2004). [Google Scholar]
• 13.Mitzenmacher M., A survey of results for deletion channels and related synchronization channels. Probab. Surv.6, 1–33 (2009). [Google Scholar]
• 14.Li R., New developments in coding against insertions and deletions. Honors thesis, Carnegie Mellon University, Pittsburgh, PA: ). [Google Scholar]
• 15.Wicker S. B., Bhargava V. K., An Introduction to Reed-Solomon Codes (Wiley-IEEE, 1994). [Google Scholar]
• 16.Vigenère B. de, Traicté des chiffres ou secrètes manières d’escrire (Abel l’Angelier, Paris: 1586), ff. 46r–49v. [Google Scholar]
• 17.Hart P. E., Nilsson N. J., Raphael B., A formal basis for the heuristic determination of minimum cost paths. IEEE Trans. Sys. Sci. Cyber.4, 100–107 (1968). [Google Scholar]
• 18.Hawkins J. A., Jones S. K., Finkelstein I. J., Press W. H., Indel-correcting DNA barcodes for high-throughput sequencing. Proc. Natl. Acad. Sci. U.S.A.115, E6217–E6226 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 19.Chen Y.-C., Liu T., Yu C.-H., Chiang T.-Y., Hwang C.-C., Effects of GC bias in next-generation-sequencing data on de novo genome assembly. PloS One8, e62856 (2013). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 20.Jain M., et al. , Nanopore sequencing and assembly of a human genome with ultra-long reads. Nat. Biotechnol.36, 338–345 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
• 21.Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, ed. 3, 2007). [Google Scholar]
• 22.Klar B., Bounds on tail probabilities of discrete distributions. Probab. Eng. Inf. Sci.14, 161–171 (2000). [Google Scholar]
• 23.Press W. H., Hawkins J. A., Jones S. K. Jr, Schaub J. M., Finkelstein I. J., HEDGES error-correcting code for DNA storage corrects indels and allows sequence constraints. NCBI BioProject.https://www.ncbi.nlm.nih.gov/bioproject/631961. Deposited 12 May 2020. [DOI] [PMC free article] [PubMed] [Google Scholar]

Supplementary Materials

Data Availability Statement

The sequenced reads used in testing are available on the Sequence Read Archive (SRA) under accession numbers SAMN14897329, SAMN14897330, SAMN14897331, SAMN14897332, SAMN14897333, SAMN14897334, and SAMN14897335 (SRA Project number PRJNA631961) (23).

Oligos used in the study are as follows:

IF538:$5^{\prime }$-TCGAAGTCAGCGTGTATTGTATG-$3^{\prime }$.

IF539:$5^{\prime }$-AACACGCTTAATCGCACTCACTA-$3^{\prime }$.

The computer code used for the generation and testing of the inner HEDGES code and outer RS code is available athttps://github.com/whpress/hedges. This paper utilized two commercial C++ source code libraries: Numerical Recipes (http://www.numerical.recipes) and the Schifra Reed-Solomon Error Correcting Code Library (http://www.schifra.com). The specific routines used in this paper are freely available for noncommercial use and are included in the above GitHub repository.