Introduction↩︎

The exponential growth of digital information has made conventional storage technologies increasingly unsustainable [1], motivating research into alternative storage media, including the use of biological molecules. DNA has emerged as a particularly promising medium owing to its exceptional density and durability [2], [3]. One gram of DNA can theoretically store exabytes of data, and DNA molecules can remain stable for centuries at room temperature under suitable conditions [2], [3], enabling compact and sustainable long-term storage. The DNA storage pipeline consists of a writing process, a storage phase, and a reading process. During writing, binary information is digitally encoded into quaternary sequences over the four DNA nucleotides Adenine (\(\mathsf{A}\)), Guanine (\(\mathsf{G}\)), Cytosine (\(\mathsf{C}\)), and Thymine (\(\mathsf{T}\)), which are then synthesized as short DNA molecules known as oligonucleotides. The encoded sequences must be short to adhere to the length limitations imposed by current DNA synthesis technologies. During reading, the stored oligonucleotides are amplified, sequenced, and the resulting reads are digitally processed and decoded to retrieve the original information.

While several proof-of-concept experiments have already demonstrated the feasibility of this pipeline, scalability remains a central challenge. The main obstacles are related to reliability, speed, and cost. DNA synthesis, amplification, and sequencing are inherently noisy biochemical processes that introduce errors and biases throughout the pipeline [4], [5]. Existing prototypes predominantly rely on high-fidelity synthesis and sequencing technologies, such as material deposition-based synthesis (Twist Bioscience) and Illumina sequencing, to mitigate these effects. Such technologies, however, are slow and expensive, limiting current systems to small-scale demonstrations that store megabytes of data rather than the exabytes envisioned. Faster and lower-cost alternatives, such as photolithographic synthesis, offer more scalable solutions but exhibit significantly higher error rates [6], [7]. A key step toward scalability is therefore to address these errors and biases algorithmically rather than preventing them biochemically through high-fidelity technologies, thereby enabling the use of faster, cheaper, yet more error-prone technologies. Achieving this goal requires the development of efficient encoders and decoders (codecs) specifically adapted to the DNA storage channel, whose distinctive noise characteristics differ fundamentally from those encountered in conventional communication and storage systems [4].

Errors in DNA storage manifest as insertion, deletion, and substitution (IDS) errors at the base level, affecting individual nucleotides [4]. In addition, biases introduced along the storage pipeline, such as those arising during PCR amplification, together with the stochastic nature of sequencing, can lead to sequence dropouts, i.e., the complete loss of encoded sequences [5]. A natural algorithmic approach to mitigate both errors and dropouts is the use of error-correcting codes (ECCs) [8], [9]. ECCs introduce structured logical redundancy into the stored data to enable reliable information retrieval under noisy conditions. A common design strategy is to employ a two-layer coding architecture in which the inner code adds redundancy within each encoded sequence to correct base-level errors, while the outer code introduces redundant sequences to compensate for dropout events [10].

ECCs provide several key benefits that enhance DNA storage systems. First, they improve data retrieval reliability by reducing the probability of decoding failures in the presence of errors and dropouts [11]. Second, they reduce the number of sequencing reads required for successful decoding, which in turn decreases both the data retrieval time and the overall sequencing effort [12]. Third, in addition to decreasing sequencing depth, ECCs also allow reducing the number of physical molecular copies required for reliable data retrieval, thereby enabling higher storage densities [13]. Fourth, they enable the correction of synthesis errors that may appear systematically across reads and therefore cannot be mitigated simply by increasing sequencing depth [14].

Due to the established benefits of ECCs, the design of codecs for DNA storage with robust error-correction capabilities has been the subject of extensive research. The ECCs adopted in existing DNA storage codecs vary widely in scope. In particular, they differ in whether they are used for error detection or error correction, whether they address base-level errors and/or dropouts, and in the types of base-level errors they can handle. Existing designs include codecs that combine dropout recovery with error detection [3], [12], [15]–[19], codecs that jointly address dropouts and substitution errors [6], [20]–[22], codecs whose ECCs support dropout recovery together with correction of all three types of IDS errors [14], [23]–[25], and codecs that focus on correcting IDS errors within individual sequences without providing explicit dropout recovery [26]. Most of these codecs also enforce content-specific constraints on the encoded DNA sequences, typically by limiting homopolymer length and balancing GC content, with some additionally avoiding specific sequence motifs [24]. Such constrained coding strategies follow the general principle of preventing or reducing errors and biases by avoiding sequence patterns believed to be problematic for the underlying biochemical processes. Some works focus exclusively on constrained coding and therefore design codecs without explicit ECCs [2], [27]–[29], including approaches that additionally account for thermodynamic constraints [27], [29]. A recent study suggests that addressing errors through ECCs is more efficient than attempting to avoid them through constrained coding alone [30].

Beyond conventional coding approaches over the standard four-nucleotide alphabet, alternative encoding strategies based on degenerate bases and composite DNA letters have also been explored [31]–[33]. Recent works have further investigated deep-learning-based reconstruction [34], neural decoding of polar codes over synchronization-error channels [35], and joint probabilistic processing tailored to Nanopore reads [36]. These approaches can exploit complex empirical error patterns and soft information from multiple reads, potentially improving robustness and reducing coding redundancy or decoding complexity. However, they typically require representative training data, accurate channel models, or specialized computational procedures, which can limit their portability across experimental workflows and complicate direct system-level comparisons.

In this work, we focus on complete DNA storage codecs based on the standard four-nucleotide alphabet and introduce a novel codec called DNA-MGC+ (Marker Guess & Check Plus). Through a broad range of in silico and in vitro settings, including varying error and bias profiles as well as experimental results obtained with both Illumina and Nanopore sequencing, we show that DNA-MGC+ consistently outperforms other selected codecs. The obtained results highlight the versatility of DNA-MGC+, demonstrating that its performance is not limited to a particular error profile or experimental workflow. Specifically, DNA-MGC+ achieves simultaneous gains across several key performance metrics under explicit reliability constraints, including the minimum required sequencing depth, read cost, decoding time, maximal error-correction capability, and storage density. These gains also persist at larger file sizes, extending beyond the practical scalability limits of the benchmarked codecs. Notable performance results include supporting reliable decoding under base-level error rates of up to 24% in synthetic scenarios, achieving an estimated storage density of 57 EB/g under experimentally derived error and bias profiles, and enabling reliable retrieval under both Illumina and Nanopore sequencing at sequencing depths below 3, with associated read costs below 3.5 nts/bit. While most previous studies evaluate codecs under high-fidelity workflows, typically involving material deposition-based synthesis (e.g., Twist Bioscience) and/or Illumina sequencing, our experimental validation also includes results combining electrochemical synthesis (GenScript) with Nanopore sequencing, both of which are substantially more error-prone than their conventional counterparts. These results show that, in addition to the performance improvements realized in typical workflows, DNA-MGC+ can also enable reliable storage in lower-fidelity and more error-prone settings, demonstrating its overall resource efficiency.

Beyond the quantitative performance gains discussed above, our experimental study reveals additional findings of broader and independent interest. In particular, we observe that comparable data retrieval performance can be achieved under both Illumina and Nanopore sequencing when using state-of-the-art Nanopore basecalling methods and efficient codecs, despite the higher error rates associated with Nanopore. We further find that codecs equipped with efficient ECCs can achieve similar retrieval performance with and without content-specific constraints on encoded sequences. Furthermore, we observe that less commonly considered thermodynamic constraints based on free energy have a more pronounced impact on data retrieval performance than conventionally imposed constraints such as limiting homopolymer length and balancing GC content.

Results↩︎

Conceptual design of the DNA-MGC+ codec↩︎

We outline below the main steps of the DNA-MGC+ codec, describing how a binary file is encoded into a collection of DNA sequences and subsequently decoded from sequencing data. This high-level overview is intended to convey the structure and main operating principles of the codec, while the full mathematical and algorithmic details are deferred to the Methods section.

1. Fragmentation: The input binary file is partitioned into non-overlapping short fragments. The fragment length is determined as a function of the desired target DNA sequence length and user-defined codec parameters.

2. Outer encoding: The collection of binary fragments is encoded using a Reed–Solomon (RS) code [37], which introduces inter-sequence redundancy by generating additional binary fragments. The RS code acts as an outer code, enabling recovery from dropouts (lost sequences) as well as correction of residual errors from other decoding stages.

3. Indexing and inner encoding: A unique binary index is appended to each fragment, after which the indexed fragments are encoded using the Marker Guess & Check Plus (MGC+) code [38], [39]. The MGC+ code acts as an inner code that introduces intra-sequence redundancy through structured encoding applied separately to each indexed fragment. This redundancy enables the correction of IDS errors within individual DNA sequences. The MGC+ encoding process includes an unconstrained binary-to-quaternary mapping (\(00 \leftrightarrow \mathsf{A}\), \(01 \leftrightarrow \mathsf{T}\), \(10 \leftrightarrow \mathsf{C}\), \(11 \leftrightarrow \mathsf{G}\)) and therefore produces coded DNA sequences without content-specific constraints.

4. Filtering: The codec optionally supports filtering of the resulting DNA sequences based on arbitrary content-specific constraints. These may include, but are not limited to, homopolymer length, GC content, sequence motifs, or thermodynamic properties. This capability is enabled by the outer RS code, which allows the generation of an excess pool of candidate sequences from which any subset can be selected to represent the file through puncturing.

5. Decoding: Decoding operates on noisy, unordered DNA sequences, which may correspond either to raw sequencing reads or to consensus sequences obtained after clustering. The sequences provided to the decoder are not required to match the reference sequence length, as insertions and deletions can be handled by the inner MGC+ code. Decoding proceeds by applying the inner MGC+ decoder, extracting fragment indices, decoding the outer RS code, and concatenating the recovered fragments to reconstruct the original file.

The main motivation behind this design is to combine strong base-level IDS error-correction capability at the inner code with sequence-level erasure and substitution correction at the outer code. The inner-code component is provided by MGC+, which offers more efficient correction of IDS errors than the inner codes typically used in existing DNA storage codecs, while also communicating decoding information to the outer code. The outer-code component is provided by the RS code, which can treat detected inner-decoding failures as sequence-level erasures (dropouts), thereby reducing the amount of redundancy required for recovery, while also correcting residual (undetected) substitution errors. Another distinguishing feature compared with existing designs is the ability to support arbitrary sequence constraints through filtering while retaining the erasure- and substitution- correction capability of RS codes. In contrast, previous filtering-based designs have typically relied on outer fountain codes, whose correction capability is largely limited to erasures. Finally, the decoding procedure is naturally parallelizable at both the inner and outer decoding stages, which improves scalability to larger files.

Overview of the evaluation framework↩︎

We evaluate the performance of the proposed DNA-MGC+ codec in comparison with selected codecs under three settings: (1) fully synthetic channel models, (2) experimentally derived in silico error and bias profiles, and (3) in vitro wet-lab experiments. In all settings, evaluation follows an end-to-end pipeline in which a binary input file is encoded into DNA sequences using the codec under consideration, processed through the associated channel model or experimental workflow, and decoded back into binary form.

In the first in silico setting, we consider a randomly generated binary file of fixed size 15 KB and evaluate the performance of DNA-MGC+ and the comparison codecs under multiple synthetic channel scenarios that capture different error and bias regimes. The following six codecs were selected for comparison based on their prominence in the literature, the availability of wet-lab validation, and outcomes of previous benchmarking efforts [13]: DNA-Aeon [24], HEDGES [14], DNA-Fountain [12], DNA-RS [6], [21], DNA-StairLoop [25], and an LDPC-based scheme [20]. Across all evaluated codecs, we consider a total of 16 configurations spanning three code-rate classes, with target rates of approximately 0.50 (low), 1.00 (medium), and 1.50 (high) bits per nucleotide (bits/nt). The design length of the encoded DNA sequences is chosen to be as close as possible to 150 nucleotides. We then increase the input file size beyond 15 KB to assess the system-level scalability of the best-performing codec configurations.

In the second in silico setting, we again consider a randomly generated binary file of size 15 KB and compare DNA-MGC+ against the same six representative codecs. In this case, however, the evaluation uses experimentally derived error and bias profiles based on the DT4DDS [5] digital twin for DNA storage. Finally, in the in vitro setting, we evaluate DNA-MGC+ using wet-lab experiments performed on a file of size approximately 24 KB. This evaluation includes a comparison with the two best-performing codec configurations identified from the preceding in silico evaluations. The corresponding experimental workflow and codec configurations are described in the dedicated sections.

Performance is assessed using three metrics, reported consistently across all evaluation settings. First, decoding reliability is measured through the minimum coverage depth required to retrieve the file with an exact match, subject to a reliability constraint that requires successful decoding in 50 out of 50 independent trials. Second, we report the read cost, defined as the ratio of the minimum required coverage depth to the code rate. Normalizing by the code rate penalizes lower-rate configurations and removes the bias that would otherwise arise when comparing only the coverage depth across codecs storing different amounts of information per nucleotide. This metric therefore enables a fair comparison between codecs operating at different code rates by quantifying the minimum number of nucleotides that must be read per stored information bit (nts/bit). Third, we report the average decoding time measured at the minimum reliable coverage depth. Decoding time excludes clustering, alignment, and consensus formation, which are handled separately to isolate the computational cost of the codec decoding algorithms themselves. All computations are performed on a machine equipped with an AMD Ryzen™ Threadripper™ PRO 7985WX CPU, 256 GB of available memory, and 64 CPU cores. The number of cores used is specified for each experiment, and a two-hour decoding time limit is imposed in all experiments. Unless otherwise stated, all results are obtained using a common preprocessing pipeline prior to decoding, in which noisy reads are clustered using CD-HIT [40], followed by multiple sequence alignment using Kalign [41] and consensus calling.

In silico performance under synthetic error and bias models↩︎

We start our performance analysis by evaluating the DNA-MGC+ codec in comparison with the other aforementioned codecs under a synthetic end-to-end channel model. This model captures key sources of variability in the DNA storage pipeline, including base-level errors, coverage bias across reference DNA sequences, and the stochastic nature of the sequencing process. The channel is parameterized by three independent quantities: a bias parameter, a coverage depth, and an error rate (Fig. 1a). Given a set of reference DNA sequences produced by the encoder, the channel generates a collection of noisy reads through three successive stages:

1. Bias: Each reference sequence is assigned a multiplicative weight drawn from a lognormal distribution with unit mean and standard deviation governed by the channel bias parameter. These weights are normalized to form a probability vector that defines the sequence-dependent sampling probabilities to be used in the next stage.

2. Coverage: Reads are generated by random sampling with replacement according to the probability vector obtained in the previous stage, where the total number of reads is equal to the number of reference sequences multiplied by the channel coverage depth parameter.

3. Error: Base-level errors are introduced randomly and independently into each read according to the channel error rate, with substitutions, deletions, and insertions occurring in fixed proportions of 53%, 45%, and 2% of the total error rate, respectively, consistent with ratios reported in the literature [5].

The coverage depth parameter fixes the average number of reads per reference sequence, while the bias parameter \(\sigma\geq0\) influences how unevenly these reads are distributed across sequences. The adopted bias model is motivated by experimental observations showing that the coverage distribution is positively skewed and is well-approximated by a lognormal distribution [5], [42]. In our simulations we consider three representative bias regimes: \(\sigma = 0\) (no bias), \(\sigma = 0.5\) (moderate bias), and \(\sigma = 1\) (strong bias). When \(\sigma = 0\), all sequences are sampled with equal probability, leading to a symmetric coverage distribution, whereas for \(\sigma > 0\) the increasing standard deviation induces progressively skewed coverage distributions (see Fig. 1b). An important consequence of skewed coverage is an increased occurrence of sequence dropouts, whereby some reference sequences receive zero reads at the channel output. The dropout rate grows with increasing bias strength (see Fig. 1c).

Performance across multiple channel scenarios↩︎

We first consider a fixed file size of 15 KB and evaluate DNA-MGC+ and the comparison codecs across 12 different synthetic channel scenarios, obtained by combining the three bias regimes with four error rates (1%, 5%, 10%, and 15%). For each scenario, we determine the minimum coverage depth required for reliable decoding for all 16 configurations of the comparison codecs, in addition to three configurations of DNA-MGC+ corresponding to target code rates of 0.50, 1.00, and 1.50 bits/nt. This is done by performing a binary search over coverage depths in the range [1, 32] with a resolution of 0.25, under the reliability constraint of 50 successful decodings (with an exact match) out of 50 trials. The corresponding read cost is then computed by normalizing the minimum coverage depth by the associated code rate. Figs. 1d–g report, for each of the 12 scenarios, the three best-performing codecs in terms of read cost, together with their corresponding coverage depths.

The results in Figs. 1d–g show that the strongest performers across the synthetic channel conditions are DNA-MGC+, HEDGES, DNA-Aeon, and DNA-RS. In particular, DNA-MGC+ achieves the lowest read cost in 11 out of the 12 bias–error combinations, with the only exception occurring under strong bias at 1% error rate, where the low-rate configuration of DNA-RS attains a lower read cost. Beyond read cost, DNA-MGC+ also achieves the lowest coverage depth among all codecs for error rates of 5%, 10%, and 15% across all bias regimes, demonstrating that its reliability advantage persists as both dropouts and errors increase. Notably, under the highest error rate of 15%, only 2 out of the 19 tested codec configurations meet the reliability constraint within the simulated coverage depth range of [1, 32]. These two cases correspond to the low-rate configurations of HEDGES and DNA-MGC+ (see Fig. 1g). Both succeed in the no-bias and moderate-bias regimes, while DNA-MGC+ is the only codec that achieves reliable decoding under strong bias at 15% error rate.

In Figs. 1h–i, we highlight results for the moderate-bias, 5% error-rate scenario in terms of average decoding time and the relation between write cost and read cost. The average decoding times reported in Fig. 1h are measured at the minimum coverage depth required by each codec for reliable decoding, with the corresponding decoding computations restricted to the use of a single CPU core for all codecs. The results show that DNA-MGC+ consistently ranks among the three fastest codecs across all three code-rate classes. Fig. 1i illustrates the trade-off between write cost and read cost for the same scenario, where the write cost is defined as the reciprocal of the code rate expressed in nts/bit. Codecs with lower write cost (that is, higher code rate) generally require higher read cost to achieve reliable decoding, reflecting the classical trade-off between logical redundancy added during encoding and redundancy supplied through coverage depth. DNA-MGC+ achieves a favorable balance in this trade-off, as it lies near the lower envelope of the cloud of points, combining relatively low write cost with some of the lowest read costs among the tested codecs.

We further examined the robustness of the low-rate configuration of DNA-MGC+, which, as shown in Fig. 1g, supports reliable decoding at error rates up to 15% even under strong bias. Building on this observation, we assessed how far the error rate can be increased while still satisfying the reliable decoding constraint within the coverage depth range of [1, 32]. Using the same preprocessing pipeline applied across all codecs (CD-HIT + Kalign), we found that reliable decoding remains achievable up to an error rate of 16% under moderate bias. Fig. 1j shows the reliable decoding region for the moderate-bias regime. Similar figures for the no-bias and strong-bias regimes are provided in Supplementary Figs. S5 and S6, showing that the maximum achievable error rate is also 16% in the no-bias case, and 15% in the strong-bias case. Our analysis further reveals that at these high error rates, clustering becomes the limiting factor as most reads end up being unclustered. This in turn renders the subsequent alignment stage ineffective since the majority of clusters contain only a single read. Adjusting the parameters of CD-HIT did not help mitigate this effect.

To address the clustering bottleneck, we tested an alternative clustering algorithm introduced by Rashtchian et al. [43], referred to here as “Clustering Billions of Reads” (CBR), combined with the same alignment method Kalign. This approach extends the reliable decoding region to 21% under moderate bias, as illustrated in Fig. 1j. For error rates of 10% or less, CBR offers no measurable advantage over CD-HIT, and the improvement becomes apparent only as the error rate approaches the regime where CD-HIT no longer clusters effectively. We also tested pairing CBR with a more computationally intensive alignment algorithm, MUSCLE [44], which further extends the reliable decoding region to 23% under moderate bias. Beyond this point, even CBR becomes ineffective as almost all reads remain unclustered. For the no-bias and strong-bias regimes, the results provided in Supplementary Figs. S5 and S6 indicate that by using CBR and MUSCLE, DNA-MGC+ achieves reliable decoding up to an error rate of 24% in the no-bias case, and 21% in the strong-bias case.

Additional results obtained under the same synthetic channel model are provided in the Supplementary Information. These include the encoding time required by each codec for the 15-KB file (Supplementary Fig. S4), results similar to those in Fig. 1 under equal insertion, deletion, and substitution proportions (Supplementary Figs. S1–S3), a comparison between DNA-MGC+ and the Gungnir codec [26] (Supplementary Note S4), and full numerical results for all 19 codec configurations across the 12 channel scenarios (Supplementary Tables S6–S8).

Scalability across file sizes↩︎

While the preceding results evaluate codec performance for a fixed 15-KB file across different channel conditions, we next assess the system-level scalability of selected codecs by varying the input file size from a few kilobytes to hundreds of megabytes. For this experiment, we focus on the moderate-bias, 5% error-rate scenario of the synthetic channel model and consider the four codecs that achieved the best performance in terms of read cost under this scenario, according to the results in Fig. 1 and Supplementary Table S7, namely DNA-MGC+ (low-rate), DNA-Aeon (medium-rate), DNA-RS (medium-rate), and HEDGES (medium-rate). The encoded sequence length is kept fixed at approximately 150 nts while the file size is varied. For DNA-MGC+, the redundancy allocation of the low-rate configuration is slightly adjusted for files larger than 50 KB by reducing the inner-code redundancy, allowing the encoded sequences to accommodate the additional packet index required for larger files while maintaining the same sequence length, as described in the Methods section and Supplementary Table S1.

Since the objective of this experiment is to evaluate system-level scalability, we allow each codec to use the full 64 CPU cores whenever supported by its implementation, unlike the preceding experiments, in which computations were intentionally limited to a single CPU core. The two-hour decoding time limit per file is nevertheless maintained. To avoid making clustering the limiting factor in this larger-scale setting, all codecs are evaluated using the CBR clustering algorithm, which scales more favorably to large numbers of reads than CD-HIT, followed by Kalign alignment and consensus calling. For each file size and codec configuration, we report the encoding time (Fig. 2a), the minimum coverage depth required for reliable decoding under the same 50-out-of-50 reliability criterion used throughout this work and the associated read cost (Fig. 2b), and the average decoding time measured at the minimum coverage depth (Fig. 2c). When a codec reaches an implementation-specific limitation or exceeds the two-hour decoding time limit, larger file sizes are not evaluated for that codec.

The results in Fig. 2 show that DNA-MGC+ scales more favorably to larger files than the comparison codecs while maintaining both reliable and efficient decoding. In particular, DNA-MGC+ successfully decodes files of up to 200 MB within the imposed two-hour decoding time limit while attaining a read cost below 5 nts/bit throughout the evaluated range. Fig. 2b further shows that the minimum coverage depth required for reliable decoding and the corresponding read cost remain nearly constant for DNA-MGC+ as the file size increases. The only exception is a small increase for file sizes above 50 KB, which is expected since the reduction in inner-code redundancy required to accommodate the additional packet index slightly reduces the overall error-correction capability. Nevertheless, DNA-MGC+ consistently achieves the lowest read cost among the evaluated codecs across all tested file sizes. Fig. 2c further indicates that the decoding time of DNA-MGC+ grows nearly linearly with file size, with a slope (growth rate) of approximately one. Both the decoding time and its growth rate are lower than those of the comparison codecs. In terms of encoding time (Fig. 2a), DNA-MGC+ exhibits performance comparable to DNA-Aeon and DNA-RS, depending on the file size, whereas HEDGES provides the fastest encoding.

The remaining codecs exhibit more limited scalability. DNA-Aeon maintains constant coverage requirements and read cost across the evaluated file sizes but reaches the two-hour decoding time limit for file sizes beyond 500 KB. DNA-RS shows a progressive increase in read cost as the file size grows from 50 KB to 500 KB, indicating that additional reading effort is required to maintain reliable decoding for larger files, before also reaching the decoding time limit. HEDGES exhibits a similar increase in read cost as the file size increases from 10 KB to 500 KB, beyond which no further experiments could be performed because of an inherent limitation rather than decoding complexity. Specifically, the available implementation of HEDGES contains hard-coded limits of 256 packets and 256 sequences per packet, which are insufficient for larger files.

In silico performance under experimentally derived error and bias profiles↩︎

Our second in silico performance evaluation follows a standardized workflow based on DT4DDS [5], a digital twin for DNA data storage built on real experimental datasets and designed to reproduce the behavior of state-of-the-art storage pipelines. In this study, we focus on the “low-fidelity” workflow defined in Ref. [13], which models a pipeline consisting of electrochemical synthesis (e.g., GenScript), amplification using an error-prone polymerase (e.g., Taq), and short paired-end Illumina sequencing. This scenario is characterized by higher error rates and more pronounced bias compared with high-fidelity pipelines based on material deposition synthesis (e.g., Twist) and low-error polymerases (e.g., Q5). As such, it provides a particularly relevant setting for assessing the ability of the codecs to compensate for the error and bias introduced by low-fidelity technologies.

Compared with the synthetic model used previously, this experimentally derived model imposes a more nuanced error profile. While the overall end-to-end error rate is fixed at approximately 1.6%, errors exhibit correlation across neighboring bases, position-dependent behavior with increased rates near sequence extremities, and asymmetric substitution patterns derived from empirical measurements. In addition, each stage of the workflow is modeled separately, such that the resulting error composition reflects the successive biochemical processes rather than a simple end-to-end mixture. For example, under the considered workflow, deletions predominantly originate during synthesis, whereas substitutions mainly arise during sequencing. More specifically, the error profile includes a 1.56% deletion rate during synthesis, a 0.011% substitution rate during amplification, and substitution rates of 0.11% and 0.25% in the forward and reverse Illumina reads, respectively, which are subsequently merged using NGmerge [45].

The experimentally derived workflow also provides a more refined treatment of coverage. Instead of a single coverage depth parameter, the model enables separate parametrization of physical redundancy and sequencing depth. Physical redundancy is defined as the average number of physical DNA molecules per reference sequence during storage, whereas sequencing depth corresponds to the average number of sequencing reads obtained per encoded sequence. This distinction is important because physical redundancy directly constrains the achievable storage density (SD), which can be approximated as \[\text{SD} \approx \frac{\text{code rate}}{\text{physical redundancy}}\times 113.7,\] where the rate of the codec is expressed in bits/nt and the storage density in exabytes per gram of DNA (EB/g) [13]. A codec that maintains reliable decoding at low physical redundancy therefore unlocks higher storage densities, since fewer molecular copies are required for data retrieval.

Using the low-fidelity workflow, we evaluate the same 19 codec configurations examined under the synthetic model on a randomly generated 15-KB file, together with an additional DNA-MGC+ configuration of code rate 0.50 bits/nt whose parametrization is optimized for this specific workflow. In the synthetic model evaluation, DNA-MGC+ was assessed using fixed parametrizations across 12 different error and bias conditions. In general, it is possible to optimize the parametrization of DNA-MGC+ according to the characteristics of the model under consideration. This corresponds to a design tradeoff between the redundancy allocated to the inner and outer codes, which we elaborate on in the Discussion section. The corresponding parameter settings for all evaluated configurations are provided in the Supplementary Information.

In one setting, we fix the physical redundancy to 100 and measure the decoding success rate as a function of sequencing depth (Fig. 3a), as well as the minimum sequencing depth and associated read cost required to satisfy the same 50-out-of-50 decoding reliability criterion used earlier (Fig. 3b). Notably, only 9 out of the 20 tested codec configurations meet the reliability constraint within the simulated sequencing depth range of [1, 32]. Among these, the low-rate, medium-rate, and optimized DNA-MGC+ configurations achieve the lowest read costs. In particular, the optimized configuration enables reliable decoding at a sequencing depth as low as 1, corresponding to a read cost of 2 nts/bit, while the two non-optimized DNA-MGC+ configurations also outperform all remaining codecs in terms of read cost. All numerical results for this setting are summarized in Supplementary Tables S9 and S10.

In another setting, we fix the sequencing depth to 30 and determine the minimum physical redundancy required for reliable decoding. While reduced physical redundancy, like reduced sequencing depth, increases the occurrence of sequence dropouts, it also introduces a more fundamental challenge: the limited number of physical molecular copies reduces variability across reads and limits the effectiveness of clustering (see Discussion section for more details). Fig. 3c reports the maximum achievable storage density (SD) for each codec configuration, obtained by applying the above SD formula using the measured minimum physical redundancy. The results show that the three DNA-MGC+ configurations achieve the highest storage densities, with the optimized DNA-MGC+ configuration enabling reliable decoding at a physical redundancy of only 1, corresponding to a storage density of approximately 57 EB/g.

In vitro experimental evaluation↩︎

To evaluate the performance of DNA-MGC+ under in vitro experimental conditions, we encoded a file of size 24 KB¹ using four different configurations of DNA-MGC+ and one configuration of each of the comparison codecs DNA-Aeon and HEDGES. DNA oligonucleotides corresponding to the six codec configurations were synthesized in a single pool via electrochemical synthesis performed by GenScript. The target design length of the encoded DNA sequences was set to 126 nts, below the maximum length of 170 imposed by GenScript. This choice allowed sufficient margin for the inclusion of 20-nt forward and reverse primer sequences, as well as a 4-nt unique tag to distinguish among the six configurations. Sequencing was performed using both Oxford Nanopore sequencing and paired-end Illumina sequencing.

The four DNA-MGC+ configurations considered in this experimental evaluation are based on two distinct inner and outer code parametrizations with different redundancy levels. The first parametrization yields an overall code rate of 1.03 bits/nt and is denoted as design \(\mathrm{A}\), while the second introduces additional redundancy and results in a lower code rate of 0.71 bits/nt, denoted as design \(\mathrm{B}\). For each design, we consider an unfiltered variant, in which no constraints are imposed on the content of the encoded DNA sequences, as well as a filtered variant. In the filtered variants, an excess number of candidate sequences is generated using the codec, and only those satisfying a set of constraints are retained. These constraints include restrictions on homopolymer length, GC content, and the avoidance of selected sequence motifs, together with an additional filtering step based on thermodynamic properties (Gibbs free energy). By design of the DNA-MGC+ codec, the filtering step does not affect the code rate, so the filtered variants of designs \(\mathrm{A}\) and \(\mathrm{B}\) retain code rates of 1.03 and 0.71 bits/nt, respectively. These constraints, however, do impact encoding and decoding complexity. Full details of the filtering criteria and parameter selection are provided in the Methods section and in Supplementary Table S1, respectively.

In addition to DNA-MGC+, one configuration of each of DNA-Aeon and HEDGES was included in the same experimental pool, with code rates of 1.0 and 0.61 bits/nt, respectively. These codecs were selected based on their performance in the in silico evaluations (Figs. 1 and 3), where both achieved reliable decoding and consistently ranked among the strongest performers behind DNA-MGC+. In particular, DNA-Aeon (1.0 bit/nt) was chosen for its favorable read cost performance at medium rate across multiple in silico scenarios (see Figs. 1d–g and Fig. 3b). HEDGES (0.61 bits/nt), on the other hand, was selected for its strong robustness to IDS errors at low rate (see Fig. 1g). Although DNA-RS also exhibited competitive in silico performance, its main advantage lies in the outer RS code, whose benefits are shared with DNA-MGC+ by design. At the inner code level, however, DNA-MGC+ provides stronger capabilities for correcting insertion and deletion errors. DNA-RS was therefore not included in the in vitro comparison, since its principal strengths are inherited by DNA-MGC+ and can be further enhanced through optimized redundancy allocation, as observed in Fig. 3.

The stored file was successfully decoded with an exact match under both Nanopore and Illumina sequencing, across all six codec configurations, albeit with quantitatively different performance outcomes. The performance comparison across multiple metrics is shown in Fig. 4. For each codec configuration, we report the minimum sequencing depth required for reliable decoding, determined via progressive read downsampling under the same 50-out-of-50 decoding reliability criterion used in the in silico evaluations. The downsampling is applied after the raw reads are processed to trim primers, filter payloads by length (\(\pm\)​10 nt of the design length), and demultiplex into per-configuration subsets using the 4-nt tags, as detailed in the Methods section. The corresponding read cost, expressed in nucleotides per information bit, is also reported. We further report the average decoding time for recovering the stored file, measured at the minimum sequencing depth required by each codec for reliable decoding. As in the previous results, decoding times exclude clustering, alignment, and consensus calling, which were performed using CD-HIT and Kalign.

For Nanopore sequencing, multiple Dorado [46] basecalling models were tested, revealing pronounced tradeoffs between basecalling accuracy and computational cost. The fast, high accuracy (hac), and super accuracy (sup) models exhibited basecalling runtimes differing by approximately one order of magnitude between successive models, with corresponding base-level error rates of 5.6  0.2%, 2.2  0.15%, and 1.4  0.1% (averaged across the six codec configurations), respectively. As expected, the results in Fig. 4a show that reduced basecalling accuracy leads to a systematic increase in the minimum sequencing depth required for reliable decoding, and consequently in the associated read cost. Across all three basecallers, all four configurations of DNA-MGC+ consistently require lower sequencing depth and read cost than both DNA-Aeon and HEDGES. In particular, under the dorado-sup basecaller, the filtered variant of design \(\mathrm{B}\) achieves the lowest minimum sequencing depth of 2.75, while the filtered variant of design \(\mathrm{A}\) attains the lowest read cost of 3.4 nts/bit. Furthermore, the filtered variants of DNA-MGC+ provide modest but consistent improvements over their unfiltered counterparts across all cases.

Under Illumina sequencing, we observed an average base-level error rate of 0.34  0.03%, which is lower than that of Nanopore sequencing, as expected. Overall, the performance trends under Illumina closely mirror those obtained with Nanopore. In particular, DNA-MGC+ outperforms both DNA-Aeon and HEDGES, with the filtered variant of design \(\mathrm{B}\) achieving the lowest minimum sequencing depth of 2.25, and the filtered variant of design \(\mathrm{A}\) attaining the lowest read cost of 2.4 nts/bit. A comprehensive performance comparison across Nanopore and Illumina platforms is provided in the Discussion section.

In terms of decoding time, Fig. 4b reveals clear differences across codecs, with consistent trends observed under both Illumina and Nanopore sequencing. When decoding is performed on a single core, HEDGES achieves the fastest decoding time, with the unfiltered variant of design \(\mathrm{B}\) of DNA-MGC+ following closely, while DNA-Aeon is approximately one order of magnitude slower. When parallelization across eight CPU cores is enabled, DNA-MGC+ achieves the fastest decoding performance overall, requiring less than 10 seconds on average for decoding. This improvement reflects the structure of the DNA-MGC+ decoder, which supports parallel processing at both the inner and outer code levels. The publicly available implementation of DNA-Aeon also supports parallel execution, and the corresponding parallel decoding results are therefore reported in Fig. 4b, whereas the available implementation of HEDGES does not support parallel decoding and is thus reported under single-core execution only. Another observation worth noting is that the unfiltered variants of DNA-MGC+ consistently decode faster than their filtered counterparts. This difference stems from the codec design, discussed in more detail in the Methods section, whereby filtered variants rely on a more complex decoding algorithm at the outer RS code level.

Discussion↩︎

The performance evaluations presented in this work demonstrate that DNA-MGC+ consistently outperforms the selected codecs across a broad range of in silico frameworks and in vitro experimental conditions. These evaluations span different error and bias profiles as well as multiple sequencing technologies with distinct basecalling models, highlighting the versatility of DNA-MGC+ in delivering superior performance across diverse operating conditions. The reliability and resource-efficiency of DNA-MGC+ are reflected in several key performance metrics. Specifically, the codec achieves simultaneous gains in coverage depth requirements, read cost, decoding time, maximal error-correction capability, and storage density. In practice, these gains translate into improved data reliability, reduced sequencing effort, and faster data retrieval. Importantly, the results further indicate that DNA-MGC+ is efficient even under low-fidelity biochemical processes. Such processes are more cost-effective, and are therefore expected to play a central role in overcoming current scalability bottlenecks in DNA data storage in the future. These improvements are a direct consequence of the coding strategy implemented in DNA-MGC+, as discussed next.

The design of the DNA-MGC+ codec follows a two-layer architecture, in which logical redundancy is introduced at both the inner and outer coding stages with complementary roles. At the inner level, an MGC+ code [39] protects against base-level errors by enabling the correction of insertions, deletions, and substitutions (IDS) within individual DNA sequences. At the outer level, a Reed-Solomon (RS) code [37] mitigates the effects of coverage bias by enabling recovery from sequence dropouts, while also correcting residual errors that remain after inner decoding. Under a fixed target code rate, this separation of roles leads to a natural redundancy allocation tradeoff: allocating more redundancy to the inner layer than to the outer layer enhances IDS error-correction capability at the expense of reduced robustness to dropout events, and vice versa. When the underlying error and bias profile is known or estimated, this tradeoff can be optimized, as reflected by the improved performance of the optimized DNA-MGC+ configuration in Fig. 3.

Several of the comparison codecs considered in this work follow the same two-layer coding architecture as DNA-MGC+, including DNA-Aeon, HEDGES, DNA-Fountain, and DNA-RS, whose respective inner and outer components are summarized in Table 1. As observed in Figs. 1 and 3, among the codecs evaluated in this work, those based on this inner-outer architecture achieve better performance than DNA-StairLoop [25] and the LDPC-based scheme [20], which follow conceptually different design principles. DNA-Fountain is a notable exception, as its reduced performance stems from limitations of its inner code, whose functionality is deliberately restricted to detection rather than correction of substitution errors. The advantage of using an RS code as an outer code is well established in coding theory [11], [47], as RS codes are optimal in recovering from sequence dropouts with minimal redundancy, while also allowing correction of residual errors after inner decoding. Consistent with this theory, DNA-RS and DNA-MGC+ consistently emerge as the strongest performers in bias-dominated regimes such as Fig. 1d and Fig. 3. In parallel, strong performance under high-error conditions is most evident for codecs whose inner codes support correction of IDS errors. This explains why DNA-MGC+, HEDGES, and DNA-Aeon outperform schemes limited to substitution-only correction in Figs. 1f–g. Overall, the gains of DNA-MGC+ with respect to existing designs arise from the combination of (i) enhanced IDS error-correction capability at the inner MGC+ layer; and (ii) the use of an outer RS code that is both optimal for mitigating sequence dropouts and capable of correcting residual inner-decoding errors.

The in vitro results in Fig. 4 further indicate how structural differences between codecs translate into distinct performance under Illumina and Nanopore sequencing. Achieving reliable decoding with limited reads is a dual challenge, as low coverage increases both sequence dropouts and the effective error rate in consensus sequences. Accordingly, the minimum sequencing depth required by a codec reflects both its outer-code robustness to dropouts and inner-code IDS correction capability. Notably, the results in Fig. 4 show that DNA-MGC+ maintains nearly identical sequencing depth requirements under Illumina and Nanopore (dorado-sup) sequencing, despite the higher base-level error rate of the latter. This stability is attributable to the strong IDS error-correction capability of the inner MGC+ code, which effectively compensates for elevated error rates and preserves overall decoding performance. In contrast, DNA-Aeon performs substantially better under Illumina and approaches the sequencing depth levels of DNA-MGC+, indicating comparable dropout robustness but weaker IDS correction capability. HEDGES exhibits similar performance across both platforms yet consistently requires higher sequencing depths than DNA-MGC+ and DNA-Aeon, suggesting that its primary limitation lies in dropout mitigation rather than in IDS correction.

Beyond the challenge of reliable decoding at low sequencing depths, the second studied in silico framework highlights a distinct and more stringent challenge: reliable decoding at low physical redundancy, which enables higher storage densities (Fig. 3c). While reduced sequencing depth limits the total number of reads, reduced physical redundancy limits the number of independent molecular copies of each reference sequence. Although both effects increase sequence dropouts, low physical redundancy introduces a more fundamental decoding challenge that cannot be resolved by simply increasing the number of reads. Specifically, when only a small number of physical copies of each DNA sequence are present, errors introduced during early stages of the storage pipeline, such as synthesis, tend to propagate across most subsequent reads corresponding to the same reference sequence. This results in reduced read variability and the emergence of systematic error patterns, which in turn limit the effectiveness of clustering and consensus calling. Consequently, a larger share of the error-correction burden must be handled by the codec itself rather than by upstream read-processing operations. The superior IDS error-correction capability of DNA-MGC+ at the inner-code level is reflected in the results shown in Fig. 3c, where reliable decoding is achieved at lower physical redundancies than comparison codecs, thereby enabling higher storage densities.

In addition to their role in dropout recovery and residual error correction, outer codes also influence the flexibility of constrained coding strategies. For this reason, despite the aforementioned advantages of RS outer codes, Fountain-based outer codes have also attracted significant attention in the literature. A key feature of Fountain codes is their ability to generate a virtually unlimited number of encoded DNA sequences. This enables constrained coding through simple filtering, where sequences that violate the desired constraints are discarded until a sufficient number of valid sequences is obtained [12]. In the proposed DNA-MGC+ codec, we introduce a design methodology that allows this filtering paradigm to be extended to outer RS codes, as detailed in the Methods section. As a result, the DNA-MGC+ codec supports filtering based on arbitrary content-specific constraints while retaining the error-correction advantages of outer RS codes.

The in vitro results in Fig. 4a show that modest but consistent performance gains are achieved when filtering is applied to DNA-MGC+. The filtering used for this experiment imposes several concurrent constraints. These include limiting the maximum homopolymer length to four, enforcing a GC content between 45% and 55%, avoiding motifs such as di- and trinucleotide repeats, and ranking candidate sequences by their predicted folding Gibbs free energy (\(\Delta G\)), with those whose \(\Delta G\) values are closest to zero retained. The limited extent of the observed gains compared to the unfiltered case can be attributed to two main factors. First, the inherent robustness of DNA-MGC+ to errors and biases allows it to compensate for content-induced effects without the need to explicitly avoid them through constrained coding. Second, the impact of factors such as long homopolymers and high GC content on modern synthesis and sequencing technologies appears to be less detrimental than suggested in earlier works, which is in line with recent independent experimental findings [30].

Despite the limited magnitude of the gains, the consistent improvement observed with filtering indicates that some of the considered sequence properties do indeed influence the retrieval process. In particular, our statistical analysis of the experimental data reveals that \(\Delta G\) is the dominant factor contributing to the observed gains. Specifically, we identify Spearman correlation coefficients of \(+0.31\) for Nanopore (and \(+0.26\) for Illumina) between the \(\Delta G\) of unfiltered DNA-MGC+ sequences and their sequencing read counts, with p-values below \(10^{-50}\) in both cases. These statistics support the hypothesis that sequences with \(\Delta G\) closer to zero (less negative) receive more reads and are therefore easier to decode. In contrast, the corresponding correlation coefficients for GC content and maximum homopolymer length are weaker, at \(-0.21\) and \(+0.08\) for Nanopore, respectively. Partial Spearman correlations further confirm that the \(\Delta G\)–coverage association is largely independent of GC content or homopolymer length, whereas the marginal association between GC content and coverage largely disappears after controlling for \(\Delta G\). These results therefore indicate that \(\Delta G\) is the dominant driver of coverage variability in our experiments. Additional details and a more comprehensive statistical analysis are provided in Supplementary Note S3.

Methods↩︎

DNA-MGC+ codec design↩︎

0.0.0.1 Encoding.

Consider an input binary file of size \(B\) bits. The DNA-MGC+ encoder outputs \(N\) reference DNA sequences that collectively represent the file, each of length \(L_{\text{ref}}\) nucleotides. The encoding is carried out through a series of structured transformations that introduce logical redundancy in order to enable reliable data retrieval under error and dropout events. This process, illustrated in Fig. 5, comprises the following sequential steps.

1. Fragmentation: The input binary file is partitioned into \(K=\lceil B/k \rceil\) non-overlapping fragments, each of length \(k\) bits.

2. Outer RS Code: The \(K\) fragments are encoded using an outer Reed-Solomon (RS) code over the finite field \(\mathbb{F}_{2^{\ell_\mathrm{out}}}\) to introduce inter-sequence redundancy, where each RS symbol consists of \(\ell_\mathrm{out}\) bits. The outer RS code generates \(c_{\text{out}}\) redundant fragments, resulting in a total of \(N = K + c_{\text{out}}\) encoded fragments and an outer code rate of \(\rho_{\text{out}}=K/N\).

3. Indexing: Each outer-encoded fragment is prepended with a unique binary index of length \(\ell_{\text{out}}\) bits, forming indexed fragments of length \(k+\ell_{\text{out}}\) bits. The index length \(\ell_{\text{out}}\) is chosen to uniquely identify any position in the outer RS codeword.

4. Inner MGC+ Code: Each indexed fragment is individually encoded using the Marker Guess & Check Plus (MGC+) code to introduce intra-sequence redundancy. The encoding proceeds in two stages, where part of the redundancy is introduced in the binary domain, and another part is added in the DNA domain after binary-to-quaternary mapping.

   1. Binary-domain redundancy: The inner MGC+ code introduces \(c_{\text{in}}+1\) parity symbols derived from an underlying RS code over \(\mathbb{F}_{2^{\ell_{\text{in}}}}\), where each RS symbol consists of \(\ell_{\text{in}}\) bits. This adds \((c_{\text{in}}+1)\ell_{\text{in}}\) redundant bits to each indexed fragment, resulting in binary sequences of total length \(k+\ell_{\text{out}}+(c_{\text{in}}+1)\ell_{\text{in}}\) bits. The first \(c_{\text{in}}\) parity symbols are referred to as guess parities, while the final parity symbol is referred to as the check parity, reflecting their distinct roles in the decoding procedure described later.

   2. DNA-domain redundancy: Each resulting binary sequence is mapped into a quaternary DNA sequence according to \(00 \mapsto \mathsf{A}\), \(01 \mapsto \mathsf{T}\), \(10 \mapsto \mathsf{C}\), and \(11 \mapsto \mathsf{G}\). Optionally, periodic 2-nt markers “\(\mathsf{AC}\)” are inserted after every \(\ell_{\text{in}}\) DNA symbols. The check parity symbol is further encoded by mapping it into a short DNA barcode of length \(\beta\), where the used barcodes are drawn from a predefined DNA codebook with minimum edit distance \(d_{\text{min}}\).

The encoder therefore outputs a total of \(N=K+c_{\text{out}}\) reference DNA sequences, where the length of each sequence is given by \[\label{eqL} L_{\text{ref}}(k)=\frac{1}{2}\bigg(\overbrace{k}^{\text{ \shortstack{payload\\fragment}}}+\overbrace{\ell_{\text{out}}}^{\text{index}}+\overbrace{c_{\text{in}}\ell_{\text{in}}}^{\text{guess-pars}}\bigg)+2\bigg(\overbrace{\frac{k+c_{\text{in}}}{2\ell_{\text{in}}}+1}^{\text{num. of markers}}\bigg)+\overbrace{\beta}^{\text{\shortstack{check-par\\barcode}}}.\tag{1}\]

In summary, the DNA-MGC+ encoder is fully specified by the parameter set (\(k,\, \ell_{\text{in}},\, \ell_{\text{out}},\, c_{\text{in}},\, \rho_{\text{out}},\, d_{\text{min}},\, \beta, \, \text{marker mode}\)), where \(k\) denotes the payload length of each fragment (in bits); \(\ell_{\text{in}}\) and \(\ell_{\text{out}}\) denote the symbol sizes (in bits) of the inner and outer codes; \(c_{\text{in}}\) is the number of guess parity symbols generated by the inner code; \(\rho_{\text{out}}=K/(K+c_{\text{out}})\) is the rate of the outer code; \(d_{\text{min}}\) and \(\beta\) denote the minimum barcode edit distance and the corresponding barcode length used for encoding the check parity; and the marker mode specifies whether periodic DNA markers are inserted or not.

The fragment payload length \(k\) is a design parameter that is selected so that the reference DNA sequence length \(L_{\text{ref}}(k)\) (see equation 1 ) matches a target design length \(L_{\text{target}}\) as closely as possible. In general, the target length cannot always be matched exactly due to structural constraints imposed by the codec. For instance, since the inner and outer codes operate over symbols of length \(\ell_{\text{in}}\) and \(\ell_{\text{out}}\) bits, respectively, \(k\) must be divisible by both \(\ell_{\text{in}}\) and \(\ell_{\text{out}}\). Moreover, the index length \(\ell_{\text{out}}\) must also be divisible by \(\ell_{\text{in}}\) since it is included in the inner code. In addition, when markers are used, the total number of inner code symbols (excluding the check parity) must be even to ensure periodicity. The underlying RS codes also impose the standard field size constraints of \(2^{\ell_\mathrm{in}}\geq (k+\ell_\mathrm{out})/\ell_{\text{in}} + c_{\text{in}} + 1\) and \(2^{\ell_\mathrm{out}}\geq K + c_{\text{out}}\). The latter constraint is equivalent to \(N \leq 2^{\ell_\mathrm{out}}\), thereby limiting the total number of encoded DNA sequences to at most \(2^{\ell_\text{out}}\) for a given file. Consequently, either \(\ell_\text{out}\) must be chosen sufficiently large to enable the generation of enough sequences within a single encoding instance, or the file must be partitioned into smaller independently encoded blocks, resulting in a multi-packet representation as discussed in the next subsection.

Taking the aforementioned constraints into account, the fragment length \(k\) is selected as follows. Given a target design length \(L_{\text{target}}\) and fixed parameters (\(\ell_{\text{in}},\, \ell_{\text{out}},\, c_{\text{in}},\, c_{\text{out}},\, d_{\text{min}},\, \beta, \, \text{marker mode}\)), the encoder determines \(k^{\star}\) that minimizes \(|L_{\text{ref}}(k) - L_{\text{target}}|\) subject to \(L_{\text{ref}}(k) \le L_{\text{target}}\). Thus, \(L_{\text{target}}\) serves as an upper bound on \(L_{\text{ref}}\), and \(k^{\star}\) minimizes the gap between the two. If an exact match is required, the remaining gap \(L_{\text{target}}-L_{\text{ref}}(k^{\star})\) is filled with random DNA symbols. The encoding steps illustrated in Fig. 5 are applied once \(k^{\star}\) is determined. If the input file size \(B\) is not divisible by \(k^{\star}\), random bits are added so that the final payload fragment has length exactly \(k^{\star}\). These padding bits are removed after decoding.

Throughout all results presented in this work, the following code parameters were fixed: \(\ell_\text{in}=8\), \(\ell_\text{out}=16\), \(d_{\text{min}}=5\), and \(\beta=12\). The different codec configurations of DNA-MGC+ were obtained by varying \(k\), \(c_{\text{in}}\), \(\rho_{\text{out}}\), and the marker mode, with the full parameter specifications for each configuration provided in Supplementary Table S1. The code rate (in bits/nt) associated with each configuration follows from the ratio \(B/(N\,L_{\text{ref}})\), where \(B\) is the input file size in bits, \(N=K+c_{\text{out}}\) is the total number of encoded DNA sequences, and \(L_{\text{ref}}\) is the length in nucleotides of each sequence.

0.0.0.2 Multi-packet representation.

For large files requiring more than \(2^{\ell_{\text{out}}}\) encoded sequences, and when increasing \(\ell_{\text{out}}\) is not practical, DNA-MGC+ uses a multi-packet representation. The input file is first partitioned into smaller parts, each of which is encoded independently using the same encoding procedure described above. The encoding workflow therefore remains the same as in Fig. 5 at the level of each packet, except that, in the indexing step, each sequence is prepended with a packet index in addition to its within-packet sequence index. The packet index is added before inner encoding and is thus protected by the inner code. This eventually allows the decoder to assign each recovered fragment to the correct packet and to reorder the decoded packets before reconstructing the original file.

If the total number of encoded sequences required for the file is \(N\), the packet size can be limited to at most \(N_P \leq 2^{\ell_{\text{out}}}\) sequences, resulting in \(P\approx N/N_P\) packets, where \(N_P\) is a user-defined codec parameter. The packet-index length is chosen to be sufficient to uniquely identify all \(P\) packets (requiring at least \(\log_2 P\) bits) while satisfying the previously described divisibility constraints imposed by the inner code. The additional packet-index bits are included in the encoded sequence length and are therefore accounted for in the computation of the fragment length \(k\) and the effective code rate. In the simulation results of Fig. 2, the packet-index length is set to \(\ell_{\text{out}}=16\) bits and the packet size is limited to at most \(N_P=2^{12}\) (\(<2^{\ell_{\text{out}}}=2^{16}\)) sequences, while all other results reported in this work use a single packet.

0.0.0.3 Filtering.

The DNA-MGC+ codec supports filtering of the encoded DNA sequences according to arbitrary content-specific constraints prior to synthesis and storage. These constraints can be specified independently of the coding structure and may include limits on homopolymer length, bounds on GC content, exclusion of specific sequence motifs, or other biochemical and thermodynamic requirements. This capability is enabled by the outer RS code, which allows the encoder to generate an excess pool of candidate sequences from which any subset can be selected to represent the file. In particular, by setting the outer redundancy parameter \(c_{\text{out}}\) to its maximal value, the encoder can generate up to \(N_{\text{max}}=2^{\ell_{\text{out}}}\) candidate sequences. These candidate sequences are then screened, and a subset of \(N\leq N_{\text{max}}\) sequences satisfying the desired constraints is retained for synthesis and storage, while the remaining sequences are discarded. If the file is large and a multi-packet representation is used, each packet is screened independently.

A limitation of this approach is that, when the constraints are strict, the number of admissible sequences among the \(N_{\text{max}}\) candidates within a given packet, denoted by \(N_{\text{adm}}\), may be smaller than the required number of encoded sequences. One possible workaround is to increase \(\ell_{\text{out}}\), thereby increasing \(N_{\text{max}}\). However, this also increases the encoded sequence length (see Fig. 5), which may make some constraints harder to satisfy. A more effective solution is to limit the size of each packet such that \(N_P \leq N_{\text{adm}}\). Since the number of packets scales as \(P\approx N/N_P\), reducing \(N_P\) increases the number of packets and the corresponding packet-indexing overhead. Nevertheless, this additional redundancy is typically modest, since the required packet-index length grows logarithmically with the number of packets.

From a coding-theoretic perspective, the filtering step corresponds to a customized puncturing of the outer RS code, since certain codeword coordinates are deliberately omitted. In the unfiltered case, outer RS decoding is performed using the standard Berlekamp-Massey [48] algorithm. When filtering is applied, decoding is carried out using the Welch–Berlekamp [49] algorithm, which is more suitable for punctured RS codes. Filtering does not affect the code rate provided that enough admissible sequences are available among the \(N_{\text{max}}\) candidates; otherwise, the effective code rate may decrease slightly due to the additional packet-indexing overhead. The main cost of filtering lies in the increased decoding complexity, which is \(O(N^{3})\) under Welch-Berlekamp, compared to \(O(N^{2})\) under Berlekamp-Massey.

0.0.0.4 Decoding.

The DNA-MGC+ decoder takes as input either noisy sequencing reads or, when applicable, consensus sequences obtained after clustering, and attempts to reconstruct the original binary file. A high-level description of the sequential decoding stages is provided below.

1. Inner MGC+ Code: The inner MGC+ decoder processes each DNA sequence independently with the objective of correcting IDS errors and recovering the indexed binary fragments prior to outer decoding. Insertions and deletions induce symbol offsets, causing misalignment of the inner code symbol boundaries. The decoder therefore seeks to estimate an offset pattern, denoted by \(\boldsymbol{\delta}\), reflecting the boundary shifts along the sequence. Given an estimate of this pattern, symbols affected by insertions or deletions are treated as erasures, while synchronized symbols are retained. In this manner, correction of IDS errors is reduced to erasure-and-substitution decoding over the underlying RS code.

   1. DNA-domain processing: The barcoded check parity is first decoded via minimum-distance decoding applied to the last \(\beta\) DNA symbols. If markers were introduced during encoding, an estimate of the offset pattern is then derived by solving a marker-centric maximum a posteriori (MAP) problem of the form \(\hat{\boldsymbol{\delta}} = \arg\max_{\boldsymbol{\delta}} P(\boldsymbol{\delta}\mid \text{marker observations})\), as detailed in Supplementary Note S1. The DNA symbols at the estimated marker locations are subsequently removed, and the remaining sequence is mapped to binary according to \(\mathsf{A}\mapsto 00\), \(\mathsf{T}\mapsto 01\), \(\mathsf{C}\mapsto 10\), and \(\mathsf{G}\mapsto 11\).

   2. Binary-domain processing: If an estimate \(\hat{\boldsymbol{\delta}}\) is available, the binary sequence is parsed accordingly. Symbols with nonzero estimated offsets are marked as erasures, while the remaining symbols are retained with their recovered values. All retained symbols, excluding the check parity, are then passed to the underlying RS decoder. The decoded fragment is subsequently validated using the check parity. If validation fails, a refined set of candidate offset patterns is generated and tested via a neighborhood search around \(\hat{\boldsymbol{\delta}}\), as described in Ref. [39]. If none of the candidates yields a valid result, the sequence is discarded. In the absence of markers, no initial estimate \(\hat{\boldsymbol{\delta}}\) is available; candidate offset patterns are therefore generated by exhaustive guessing and tested in turn, resulting in increased decoding complexity, as discussed in Refs. [38], [50].

2. Outer RS decoding: Indices for which no valid fragment is recovered at the inner decoding stage result in erasures at the outer-code level, while incorrectly decoded fragments induce residual substitution errors that remain undetected. All recovered fragments are reordered according to their prepended indices, and outer RS decoding is then applied. Assuming the data is encoded in a single packet, let \(E\) and \(S\) denote the numbers of erasures and substitution errors, respectively, among the \(N\) outer code symbols. Successful data retrieval with an exact match is guaranteed whenever \(E + 2S \le N - K\), which is equivalent to requiring that the number of correctly decoded fragments exceeds the number of incorrectly decoded fragments by at least \(K\). In the multi-packet case, the packet index is additionally used to first assign recovered fragments to their respective packets before outer RS decoding is applied independently within each packet. Consequently, successful data retrieval requires that the aforementioned condition be satisfied for every packet.

3. Defragmentation: The recovered payload fragments are concatenated to reconstruct the original binary file.

Computational environment and software↩︎

All simulations and decoding experiments were conducted on a machine equipped with an AMD Ryzen™ Threadripper™ PRO 7985WX CPU and 256 GB RAM. The minimum required coverage (or sequencing) depth was determined via binary search over the interval [1, 32] with a resolution of 0.25, under the reliability constraint of 50-out-of-50 successful decodings with an exact file match. Reported decoding times correspond to the wall-clock time measured at the minimum coverage depth. Clustering, alignment, and consensus calling were excluded from decoding time measurements. All reported decoding times correspond to using a single CPU core, except for Fig. 4b, where results using 8 parallel CPU cores are also included.

Read clustering was performed using CD-HIT [40] (v4.8.1) or CBR [43], and multiple sequence alignment was carried out using Kalign [41] (v3.4.0) or MUSCLE [44] (v5.1), as specified for each experiment. For CD-HIT, the identity threshold and word size were set to 0.80 and 5, respectively, with all other parameters left at their default values. For CBR, we used the implementation provided as part of the DNAStorageToolkit [51], modifying only the clustering diameter parameter to 50. To simulate the different codec configurations of DNA-Aeon, HEDGES, DNA-Fountain, and DNA-RS, we used the dt4dds-benchmark [13], [52] (v1.0.0) benchmarking suite. For DNA-Stairloop [25] and the LDPC-based scheme [20], we used the respective GitHub repositories referenced in those works. For DNA-MGC+, we used our Python implementation, referenced in the Code Availability statement, requiring Python 3.9 or later and relying only on numpy (v1.21), scipy (v1.7), reedsolo (v1.5), and galois (v0.0.3).

The channel model used in the first in silico setting was implemented as part of this work and is included in the provided code. The second in silico setting is based on the DT4DDS [5] digital twin framework and corresponds to the low-fidelity workflow described in Ref. [13]. We used the implementation of this workflow provided in dt4dds-benchmark [52] (v1.0.0) as a black-box simulator to generate sequencing reads under specified physical redundancy and sequencing depth. For the wet-lab experiment, the primer sequences were generated using the DSP Tools software [53]. The folding Gibbs free energy (\(\Delta G\)) was computed using the NUPACK Python API[54], [55] (v4.0.1.8) with default parameters at \(55^{\circ}\mathrm{C}\). Nanopore basecalling was performed with Dorado [46] (v1.3.0, simplex mode) using the R10.4.1 E8.2 400 bps fast, hac, and sup models (v4.3.0), with a minimum Q-score threshold set to 7. Paired-end Illumina reads were merged using FLASH [56] (v1.2.11). Primer trimming and payload extraction from raw sequencing reads were performed using Cutadapt [57] (v5.2), with the maximum mismatch rate set to 15%. The base-level error rates reported for the in vitro results were computed by aligning the reads to their corresponding reference sequences using minimap2 [58] (v2.30).

In vitro experimental workflow↩︎

0.0.0.5 Stored file and codec configurations.

A compressed 24 KB text file (zlib-compressed Universal Declaration of Human Rights in four languages) was encoded using four configurations of DNA-MGC+ and one configuration each of DNA-Aeon and HEDGES. For DNA-MGC+, two parameter sets were considered, referred to as Design \(\mathrm{A}\) and Design \(\mathrm{B}\). For each design, both unfiltered and filtered variants were evaluated. The corresponding number of reference sequences, payload length, and code rate are given in Table 2. Full parameter values for all evaluated codecs are provided in Supplementary Tables S1–3.

0.0.0.6 Oligonucleotide design and synthesis.

For each codec configuration, the encoded payload sequences were right-padded with random nucleotides to obtain a common length of 126 nts. A configuration-specific 4-nt tag was then appended to the right of each sequence (see Table 2), followed by the addition of a 20-nt forward primer \(\mathsf{AGCGTGCGTTACTTAGATAC}\) and a 20-nt reverse primer \(\mathsf{TCACCGTATTGCGTAGTATG}\). The resulting sequences had a total length of 170 nt, matching the maximum length supported by the synthesis provider GenScript. A single oligo pool containing all six configurations was synthesized.

0.0.0.7 Constraint-based filtering.

For the filtered DNA-MGC+ variants, an excess pool of candidate sequences was first generated using the codec. The encoded payloads were then formatted to match the 170-nt oligonucleotide structure described above, including the configuration-specific tags and primers. The sequences were subsequently screened and retained only if they satisfied the following constraints: maximum homopolymer length \(\leq 4\), GC content between 45% and 55%, and absence of di- and trinucleotide repeat motifs with four or more repetitions. For the chosen value of \(\ell_{\text{out}}=16\), the DNA-MGC+ codec allows the generation of up to \(2^{\ell_{\text{out}}}=65536\) candidate sequences, of which \(35815\) (\(54.65\%\)) and \(42349\) (\(64.62\%\)) satisfied all three constraints for designs A and B, respectively. For the retained sequences, the folding Gibbs free energy (\(\Delta G\), in kcal/mol) was computed at \(55^{\circ}\mathrm{C}\), corresponding to the primer annealing temperature (\(T_a\)). The \(1532\) and \(2277\) sequences for designs A and B, respectively, whose \(\Delta G\) values were closest to zero were then preferentially selected. Histograms of the maximum homopolymer length, GC content, and \(\Delta G\) for the encoded sequences of each of the six codec configurations are shown in Supplementary Fig. S7.

0.0.0.8 PCR, library preparation, and sequencing.

The oligo pool from GenScript, received dry, was resuspended to a concentration of 4.5 ng \(\mu\)L\(^{-1}\) using nuclease-free water. For PCR amplification, 2 \(\mu\)L of the resuspended oligo pool (4.5 ng \(\mu\)L\(^{-1}\)) was mixed with 25 \(\mu\)L KAPA HiFi HotStart ReadyMix 2\(\times\) (Roche), 1 \(\mu\)L of 10 \(\mu\)M forward primer, 1 \(\mu\)L of 10 \(\mu\)M reverse primer, and 21 \(\mu\)L nuclease-free water, yielding a final reaction volume of 50 \(\mu\)L. The reaction was prepared in three separate tubes. Thermocycling followed established protocols, consisting of an initial denaturation at \(95^{\circ}\)C for 5 min, followed by 30 cycles of denaturation at \(95^{\circ}\)C for 20 s, annealing at \(55^{\circ}\)C for 30 s, and extension at \(72^{\circ}\)C for 20 s. PCR products were subsequently purified using SPRI magnetic beads (Beckman Coulter) and analyzed on a BioAnalyzer system using the High Sensitivity DNA kit and reagents (Agilent Technologies).

Libraries for Oxford Nanopore sequencing were prepared using the SQK-LSK114 ligation sequencing kit (Oxford Nanopore Technologies). Nanopore sequencing was performed on a PromethION platform using R10.4.1 PromethION flow cells. Libraries for Illumina sequencing were prepared using the NEBNext UltraExpress^® DNA Library Prep Kit (New England Biolabs). Illumina sequencing was performed on a NextSeq 2000 system (Illumina) using a P2 flow cell and a 600-cycle kit to generate paired-end \(2\times300\) bp reads. Coverage distributions derived from the Illumina and Nanopore sequencing data for each of the six codec configurations are shown in Supplementary Fig. S8.

0.0.0.9 Read processing and decoding.

Raw sequencing reads were first processed to locate and trim the primer sequences and to extract the payload in between. For Illumina data, paired-end reads were merged prior to primer trimming, whereas Nanopore reads were processed directly. Reads with an extracted payload length within 10 nucleotides of the design length of 130 nts, that is between 120 and 140 nts, were retained for further analysis. The retained payload sequences were then demultiplexed into six separate files according to their configuration-specific 4-nt tags listed in Table 2. Specifically, a read was assigned to a codec configuration if the last four nucleotides of its extracted payload exactly matched the corresponding tag.

For each demultiplexed set, clustering, alignment, and consensus calling were performed. The resulting consensus sequences were directly provided to the corresponding decoder of each codec without additional processing. To determine the minimum sequencing depth required for reliable decoding, progressive downsampling was applied to the demultiplexed reads prior to clustering. For each codec configuration and each tested depth value, a number of reads equal to the product of the depth and the corresponding number of reference sequences was randomly sampled. Clustering, alignment, consensus calling, and decoding were then repeated at each depth value until the minimum depth enabling reliable decoding was identified.

Data availability↩︎

The reference sequences and the corresponding FASTQ sequencing data generated in the in vitro experiments are available on Zenodo at https://doi.org/10.5281/zenodo.18848901.

Code availability↩︎

The source code of DNA-MGC+ is publicly available at https://github.com/ramy-khabbaz/MGCP.

Acknowledgements↩︎

This work was supported by the France 2030 investment plan, managed by the National Research Agency (ANR), through the PEPR MolecularXiv project (ANR-22-PEXM-003) and the Initiative of Excellence Université Côte d’Azur (ANR-15-IDEX-01). The authors thank the Institute of Molecular and Cellular Pharmacology (IPMC), Sophia Antipolis, France, and its director Pascal Barbry for providing access to laboratory facilities and technical support that made the experimental work possible.

Author contributions↩︎

R.K. and S.K. designed the DNA-MGC+ codec and performed simulations and data analysis. R.K. developed the implementation code and prepared the figures and tables. S.K. wrote the manuscript, with contributions from R.K. and J.M. All authors participated in the design of the wet-lab experiment. J.M. performed the experiments and devised the thermodynamic sequence filtering procedure. M.A. and S.K. supervised the study. All authors read, revised, and approved the final manuscript.

Competing interests↩︎

The authors declare no competing interests.

Supplementary Notes↩︎

Supplementary Note S1: Detailed marker-based offset estimation for inner MGC+ code↩︎

Supplementary Note S1 extends the offset estimation framework originally introduced in Ref. [39], where the MGC+ code (the inner code of the DNA-MGC+ codec) was developed and analyzed in the binary domain. While the overall marker-based synchronization and trellis-based decoding principles remain conceptually the same, the analysis below reformulates the system in the DNA domain. The key differences arise from (i) base-level insertion, deletion, and substitution (IDS) errors in the underlying channel model, which modify the transition probability structure, and (ii) the explicit derivation of received DNA marker transition probabilities under quaternary alphabet statistics (Table 3).

Drift and Trellis Formulation

For marker-based encoding, a fixed DNA marker (\(\mathsf{AC}\)) is inserted periodically after every \(\ell_{\text{in}}\) DNA symbols. At the decoder, synchronization is recovered by constructing a trellis whose states represent the cumulative synchronization drift after each marked block. Let the trellis state space be defined as \(\Omega = \{-\Phi, -\Phi+1, \dots, \Phi\},\) where \(\Phi\) is an upper bound on the cumulative drift induced by insertion and deletion events. Let \(\mathbf{z} = (z_0, z_1, \dots, z_v)\) denote the drift sequence, where \(z_i \in \Omega\) is the cumulative offset after the \(i\)-th block and \(z_0 = 0\). The drift evolves according to \(z_i = z_{i-1} + d_i,\) where \(d_i\) is the offset increment associated with block \(i\). The received DNA sequence \(\mathbf{y}\) is segmented into blocks according to the hypothesized drift trajectory as \[\tilde{\mathbf{b}}^{(i)} \triangleq \mathbf{y}_{[(i-1)\ell' + z_{i-1} + 1,\; i\ell' + z_i]}, \quad i \in [v], \label{eq:block95segmentation}\tag{2}\] where \(\ell'=\ell+2\). The offset increment is equivalently expressed as \(d_i = z_i - z_{i-1}.\) The last \(\mu = 2\) nucleotides of each block correspond to the received marker sequence, denoted by \(\tilde{\mathbf{m}}^{(i)}\).

Marker-centric MAP Drift Estimation

Assuming independent and identically distributed IDS errors across nucleotide positions, the joint probability of observing the received sequence and a particular drift trajectory factorizes as \[\Pr(\tilde{\mathbf{y}}, \mathbf{z}) = \prod_{i=1}^{v} \Pr(\tilde{\mathbf{b}}^{(i)}, z_i \mid z_{i-1}), \label{eq:joint95probability}\tag{3}\] with \(z_0 = 0\). To reduce decoding complexity, drift estimation is performed using only the marker observations. This leads to the following marker-centric maximum a posteriori (MAP) objective: \[\hat{\mathbf{z}} = \arg\max_{\mathbf{z} \in \mathcal{Z}} \prod_{i=1}^{v} \Pr(\tilde{\mathbf{m}}^{(i)}, z_i \mid z_{i-1}), \label{eq:map95estimate}\tag{4}\] where \(\mathcal{Z}\) denotes the set of all admissible trellis paths constrained to the state space \(\Omega\). A graphical illustration of the resulting trellis and its admissible transitions is provided in Fig. 6.

Transition Probabilities

A trellis transition from state \(z_{i-1}\) to \(z_i = z_{i-1} + d\) is weighted by the probability of observing a received DNA marker \(\tilde{\mathbf{m}}\) together with a net offset \(d\) over the corresponding block. Since each block consists of a data segment followed by a fixed DNA marker, this probability decomposes as \[\Pr(\tilde{\mathbf{m}}, d) = \sum_{t=-\mu}^{\mu} \Pr(\mathrm{D}_{\mathrm{data}} = d - t)\, \Pr(\tilde{\mathbf{m}}, \mathrm{D}_{\mathrm{mar}} = t), \label{eq:transition95probability}\tag{5}\] where \(\mathrm{D}_{\mathrm{data}}\) and \(\mathrm{D}_{\mathrm{mar}}\) are random variables representing the cumulative offsets introduced by insertion and deletion events in the data portion and the marker portion of the block, respectively.

The data portion of each block consists of \(\ell_{\text{in}}\) DNA nucleotides. Under the standard base-level IDS channel with deletion probability \(P_d\), insertion probability \(P_i\), substitution probability \(P_s\), and \(P_r = 1 - P_d - P_i - P_s\), the offset distribution \(\Pr(\mathrm{D}_{\mathrm{data}} = t)\) depends only on the block length and channel parameters. Mathematically, it is given by \[\Pr(\mathrm{D}_{\mathrm{data}} = t) = \sum_{j=\max\{0,-t\}}^{\left\lfloor\frac{\ell_{\text{in}} - t}{2}\right\rfloor} \binom{\ell_{\text{in}}}{j,\, j+t,\, \ell_{\text{in}} - 2j - t} \times P_d^{\,j}\, P_i^{\,j+t}\, P_r^{\,\ell_{\text{in}} - 2j - t}. \label{eq:data95offset95pmf}\tag{6}\]

The marker contribution \(\Pr(\tilde{\mathbf{m}}, \mathrm{D}_{\mathrm{mar}} = t)\) depends explicitly on the known transmitted DNA marker \(\mathbf{m} = \mathsf{AC},\) the channel parameters, and the received marker sequence \(\tilde{\mathbf{m}} \in \{\mathsf{A},\mathsf{C},\mathsf{G},\mathsf{T}\}^{\ast}\). Since the marker length is short (\(\mu = 2\)), this probability can be evaluated exactly by enumerating all possible insertion, deletion, and substitution events affecting the marker nucleotides. For completeness, Table 3 reports \(\Pr(\tilde{\mathbf{m}}, \mathrm{D}_{\mathrm{mar}} = t)\) for all possible received marker sequences \(\tilde{\mathbf{m}}\) and offsets \(|t| \leq 2\), expressed in closed form as functions of \(P_d\), \(P_i\), \(P_s\), \(P_r\) and \(P_0\) where \(P_0 = 0.25\).

Dynamic Programming Solution

The MAP problem in 4 is solved efficiently using dynamic programming. Let \(\alpha_i(z_i)\) denote the maximum path metric ending at drift state \(z_i\) after processing \(i\) blocks. The recursion is \[\alpha_i(z_i) = \max_{z_{i-1} \in \Omega} \alpha_{i-1}(z_{i-1})\, \Pr(\tilde{\mathbf{m}}^{(i)}, z_i \mid z_{i-1}), \label{eq:dp95metric}\tag{7}\] with initialization \(\alpha_0(0) = 1\) and \(\alpha_0(z) = 0\) for \(z \neq 0\). The maximizing predecessor state is stored as \[\mathrm{pred}_i(z_i) = \arg\max_{z_{i-1} \in \Omega} \alpha_{i-1}(z_{i-1})\, \Pr(\tilde{\mathbf{m}}^{(i)}, z_i \mid z_{i-1}). \label{eq:dp95predecessor}\tag{8}\]

After processing all \(v\) blocks, the estimated drift sequence \(\hat{\mathbf{z}}\) is recovered by traceback starting from \[\hat{z}_v = \arg\max_{z_v \in \Omega} \alpha_v(z_v), \label{eq:final95drift}\tag{9}\] and the corresponding offset pattern is obtained as \(\hat{d}_i = \hat{z}_i - \hat{z}_{i-1},\) which is subsequently used for inner decoding.

Supplementary Note S2: Additional details and comments about some of the comparison codecs↩︎

LDPC Codec

The LDPC-based scheme [20] uses a conceptually different design from the inner–outer architectures adopted by most existing DNA storage codecs. The entire data stream is encoded using a single long binary LDPC code, whose encoded bits are then segmented into short fragments and mapped to DNA symbols. Each fragment carries a BCH-protected index and a marker to assist in handling insertions and deletions during decoding.

In our in silico evaluations, we used the implementation provided at https://github.com/shubhamchandak94/LDPC_DNA_Storage, keeping its default parameters and setting the oligo_length parameter to \(150\). To obtain different code rates, we replaced the underlying parity-check matrices using three matrices from https://github.com/shubhamchandak94/LDPC_DNA_storage_data/tree/master/matrices, resulting in three codec configurations with code rates of \(0.58\), \(0.67\), and \(0.80\) bits/nt. Among these, the configuration with rate \(0.58\) achieved the best overall performance under the synthetic channel model, as shown in the subsequent Supplementary Tables. However, none of the three configurations was able to decode successfully under the conditions used in the second in silico setting.

DNA-Stairloop

The DNA-StairLoop codec [25] employs a concatenated coding structure consisting of a row code and a column code, but differs from the classical inner–outer architecture through the use of a staircase interleaver that couples the two component codes across adjacent blocks. Decoding follows an iterative soft-information exchange between the row and column decoders. The decoder is designed to operate directly on unclustered sequencing reads and to exploit the presence of multiple reads within its own probabilistic decoding algorithm, without relying on consensus calling. To remain consistent with the decoding procedure proposed by the authors, clustering and alignment were therefore not performed for this codec in our evaluations.

In our in silico experiments, we used the reference implementation available at https://github.com/Guanjinqu/StairLoop. The parameter msg_length was fixed to \(19\) to obtain encoded DNA sequences of length \(147\), close to our desired target length of \(150\). Two values of the parameter block_num, \(34\) and \(68\), were used to vary the total number of encoded sequences, resulting in two codec configurations with code rates of \(0.98\) and \(0.49\) bits/nt, respectively. Our results indicate that this codec is highly sensitive to sequence dropouts under the considered parameters, leading to degraded performance under strong bias conditions. The configuration with rate 0.98 bits/nt outperformed the configuration with rate 0.49 bits/nt, suggesting that increasing the number of encoded sequences while keeping the sequence length fixed does not necessarily improve the performance of this codec. This observation is also implicit in Ref. [25], where the authors evaluate their codec using sequence lengths closer to 200, whereas comparison codecs are evaluated at lengths closer to 150, suggesting that DNA-Stairloop may benefit from allocating more redundancy within individual sequences, and hence operating at longer sequence lengths.

DNA-Aeon

As reflected in the results in Supplementary Tables S7–S9, the decoding time of DNA-Aeon is significantly higher than that of the other codecs evaluated in this work. Decoding also becomes substantially slower at lower code rates and higher error rates. When the medium-rate configuration of DNA-Aeon was tested at an error rate of 15%, decoding ran for more than two hours without producing a result and was therefore terminated; consequently, no results are reported for this setting. The same occurred for the low-rate configuration at error rates of 10% and 15%.

Supplementary Note S3: Discussion about the effect of \(\Delta G\) on the coverage distribution↩︎

As noted in the Discussion section of the main text, the filtered variants of DNA-MGC+ provide modest but consistent gains in terms of sequencing depth (and read cost) compared to their unfiltered counterparts. This is reflected in the coverage distributions shown in Supplementary Fig. S8, where the unfiltered variants consistently exhibit a substantially larger fraction of underrepresented reference sequences than the filtered variants. For example, under Illumina sequencing for Design B of DNA-MGC+, approximately 200 reference sequences in the unfiltered variant receive between 0 and 10% of the mean coverage, whereas this number decreases to roughly 50 sequences in the filtered variant. The same qualitative behavior is observed across both designs and under both sequencing platforms, demonstrating that filtering systematically reduces the number of poorly covered sequences. These strongly underrepresented sequences contribute to increased sequence dropouts at low sequencing depths, which explains the observed difference in the minimum sequencing depth required for reliable decoding between the filtered and unfiltered variants.

Our analysis of the per-sequence coverage shows that the decoding performance gains of the filtered variants are primarily attributable to the filtering criterion based on the Gibbs free energy \(\Delta G\). We support this claim through three complementary statistical analyses of the per-sequence coverage distribution obtained for the unfiltered DNA-MGC+ variants. The numerical results provided below correspond to Nanopore sequencing; qualitatively similar trends were observed under Illumina sequencing.

Characteristics of underrepresented sequences.

We first identify the subset of reference sequences whose read count falls below 10% of the mean pool coverage. Across the complete unfiltered pool, the mean \(\pm\) standard deviation (SD) of the maximum homopolymer length, GC content, and \(\Delta G\) are \(4.2 \pm 1.1\), \(48.7\% \pm 3.3\%\), and \(-4.1 \pm 2.3\) kcal/mol, respectively. Within the underrepresented subset, the corresponding means are \(3.9\), \(49.9\%\), and \(-6.8\) kcal/mol, yielding standardized differences from the pool-wide means of \(-0.27\) SD, \(+0.36\) SD, and \(-1.17\) SD. The mean \(\Delta G\) of the underrepresented sequences is thus shifted by approximately \(1.2\) pool-wide standard deviations toward more negative values, which is signficant on its own and also substantially larger than the shifts observed for GC content or homopolymer length. Notably, the homopolymer-length shift is small and negative, whereas a larger and positive shift would be expected if long homopolymers were a major contributor to lower coverage.

Marginal Spearman correlations.

Spearman correlation analysis further supports the dominant role of \(\Delta G\). Across the full unfiltered pool, \(\Delta G\) shows the strongest marginal association with per-sequence read count (\(\rho = +0.315\), \(p\)-value \(< 10^{-50}\)), compared with GC content (\(\rho = -0.205\)) and maximum homopolymer length (\(\rho = +0.082\)). The positive sign of the \(\Delta G\)–coverage correlation confirms that sequences with \(\Delta G\) closer to zero (less negative) tend to receive more reads. The weak and positive marginal correlation between maximum homopolymer length and coverage is inconsistent with the hypothesis that long homopolymers cause sequence underrepresentation.

Partial Spearman correlations.

To verify that the association between \(\Delta G\) and coverage is not merely a consequence of its relationship with GC content or homopolymer length, we compute partial Spearman correlations that isolate the \(\Delta G\)–coverage association after mathematically removing the contribution explained by each of the other variables. The partial correlation between \(\Delta G\) and coverage remains largely unchanged when controlling for GC content (\(\rho = +0.231\)) or for homopolymer length (\(\rho = +0.295\)), indicating that the association of \(\Delta G\) with coverage is largely independent of these two variables. By contrast, the marginal association between GC content and coverage (\(\rho = -0.205\)) decreases markedly after controlling for \(\Delta G\) (\(\rho = -0.062\)), suggesting that a substantial portion of the apparent GC–coverage association can be attributed to the correlation between GC content and \(\Delta G\) rather than to an independent effect of GC content on coverage.

Summary.

Taken together, these three statistical results indicate that \(\Delta G\) is the dominant sequence-level factor associated with coverage variability in our experiment, and that GC content and homopolymer length have little to no independent effect on sequence underrepresentation. A plausible mechanistic explanation is that sequences with more negative \(\Delta G\) tend to form stable secondary structures, which may interfere with PCR amplification efficiency and thereby reduce their representation in the sequenced pool. This is also consistent with the observation that the fraction of underrepresented sequences is more pronounced under Illumina sequencing, which involves additional PCR steps during library preparation and cluster generation compared with Nanopore sequencing (see Supplementary Fig. S8).

Supplementary Note S4: Comparison with the Gungnir codec↩︎

We provide here a comparison between DNA-MGC+ and the recently introduced Gungnir codec [26] under the synthetic error and bias models considered in the main text. Gungnir provides strong IDS error-correction capability by using a complex decoding procedure, in which candidate reconstructions are exhaustively tested until a hash signature associated with the encoded fragment is satisfied. This mechanism gives Gungnir high error tolerance; however, the codec also has two important limitations.

First, as explicitly acknowledged by the authors in Ref. [26], Gungnir does not provide robustness to sequence dropouts. As a result, successful file retrieval requires that each encoded sequence be read at least once. Second, Gungnir relies on an exhaustive hypothesis-search procedure whose computational cost becomes substantial at high error rates. For example, the authors report that decoding a 19.4-KB file required more than 37 hours at high error rate on a dual 32-core Xeon Platinum 8369C server. Therefore, to keep the comparison computationally tractable, we evaluate Gungnir on a small, randomly generated 5-KB file at a total error rate of 5%, for the same bias regimes (no bias, moderate bias, strong bias) considered in Fig. 1 of the main text. Unlike the experiments reported in Fig. 1, which restrict each codec to a single CPU core, the Gungnir decoder is allowed to use the full 64-core CPU of our workstation, equipped with an AMD Ryzen Threadripper PRO 7985WX processor. The same computing resources are provided to DNA-MGC+ for a fair comparison.

Following the evaluation framework of Fig. 1, we consider three code-rate classes. For Gungnir, the high- and medium-rate configurations are obtained using the Gungnir-Trit mode of the public implementation (https://github.com/HKU-BAL/Gungnir), with rates of 1.50 and 1.00 bits/nt, respectively. The low-rate configuration is obtained using the standard Gungnir mode, with a rate of 0.50 bits/nt. In all three cases, the encoded sequence length is set to 150 nts. For each configuration and bias regime, we determine the minimum coverage depth and associated read cost required for successful file retrieval at the 5% error rate. To reduce the computational burden, the reliability requirement is relaxed relative to Fig. 1: instead of requiring successful retrieval in 50 out of 50 independent trials, we require exact file recovery in a single trial at each evaluated coverage depth. The resulting coverage depth, read cost, and decoding time are reported in Fig. B below.

The results show that DNA-MGC+ requires lower coverage depth and lower read cost than Gungnir across the evaluated bias regimes and code-rate classes, while also achieving orders-of-magnitude faster decoding. In the strong-bias case, Gungnir does not achieve successful retrieval within the evaluated coverage-depth range of [1, 32], which is expected since the codec is not robust to dropouts. These results therefore show that, although Gungnir provides strong IDS error correction, DNA-MGC+ outperforms it under the evaluation framework considered in this work.

Supplementary Tables↩︎

6pt

6pt

6pt

6pt

6pt

6pt

Supplementary Figures↩︎

References↩︎

---

1. The file is obtained by zlib compression of a text document containing the Universal Declaration of Human Rights in four different languages.↩︎