De Bruijn Graph Partitioning for Scalable and Accurate DNA Storage Processing

De Bruijn Graph Partitioning for Scalable and Accurate DNA Storage Processing | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results De Bruijn Graph Partitioning for Scalable and Accurate DNA Storage Processing Florestan De Moor , Olivier Boullé , View ORCID Profile Dominique Lavenier doi: https://doi.org/10.1101/2025.05.19.654814 Florestan De Moor 1 IRISA-CNRS, Univ. Rennes , Inria, Rennes, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Olivier Boullé 1 IRISA-CNRS, Univ. Rennes , Inria, Rennes, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Dominique Lavenier 1 IRISA-CNRS, Univ. Rennes , Inria, Rennes, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Dominique Lavenier For correspondence: lavenier{at}irisa.fr Abstract Full Text Info/History Metrics Data/Code Preview PDF A bstract Motivation DNA-based data storage offers a compelling solution for long-term, high-density archiving. In this framework, accurately reconstructing high-quality encoded sequences after sequencing is critical, as it has a direct impact on the design of error-correcting codes optimized for DNA storage. Furthermore, efficient and scalable processing is essential to manage the large volumes of data expected in such applications. Results We introduce a novel method based on de Bruijn graph partitioning, enabling fast and accurate processing of sequencing data regardless of the underlying sequencing technology and without requiring prior knowledge of the information encoded in the oligonucleotides. Evaluated on both synthetic and real datasets, the method achieves excellent precision and recall. It is implemented in C++ within the software ConCluD and optimized for multi-core servers. Our experiments demonstrated that a dataset of 55 million reads, corresponding to a 135 MB binary file, can be processed in less than 10 minutes on a 16 hyper-threaded core server. 1. Introduction The world is experiencing an exponential surge in digital information, posing significant challenges for archival storage. Existing technologies consume large amounts of energy, occupy substantial physical space, incur high infrastructure costs, and require extensive maintenance to ensure long-term data preservation. Hence, new storage solutions need to be found. Advances in molecular synthesis and the rise of efficient sequencing technologies have opened the possibility of storing data in DNA molecules. This approach is gaining increasing interest as a potential solution to overcome the current technological barriers hindering its widespread adoption [ 1 ] [ 2 ] [ 3 ] [ 7 ]. DNA-based data storage follows several key steps. First, binary files are encoded into a DNA sequence representation. Then, DNA molecules are synthesized based on this encoded data and stored in a protected environment. To retrieve the information, the process begins by selecting the DNA molecules that correspond to the desired files. These molecules are then sequenced, and the resulting data are transmitted to the decoding stage to reconstruct the original files. Current synthesis processes can only produce short oligonucleotides with a maximum length of about 300 nucleotides. As a result, files must be divided into smaller fragments, each of which is indexed to allow the reconstruction of the original information. For example, a video file would be fragmented into hundreds of millions of oligonucleotides. In a single storage unit, thousands of different documents can co-exist. To identify them, each file can be assigned a unique nucleotide sequence as a tag. This tag not only helps distinguish files, but can also be used to biochemically select the oligonucleotides associated with a specific file. In this case, the tag can be designed as a pair of primers, which enables the targeted amplification of the molecules using PCR. In this framework, reading a file involves selecting the corresponding molecules by performing a PCR using the specific pair of primers associated with it, followed by sequencing the PCR product. However, the sequencing output cannot be sent directly to the decoding stage. The data must first be processed to eliminate redundancy introduced during sequencing and to reconstruct the initial set of sequences that represent the file as accurately as possible. The first processing step involves clustering the reads to group similar sequences together. In the second step, a consensus sequence is established within each group, allowing for the correction of potential errors. This process ensures the production of high-quality sequences, which can then be transmitted to the decoder stage. Since sequencing is a stochastic process, it is highly likely that not all sequences will be successfully reconstructed, and some of those provided to the decoder may still contain errors. To address this issue, the encoding process incorporates error-correcting codes. These codes enable the reconstruction of the complete information from only a subset of recovered sequences while also correcting residual errors, ensuring data integrity [ 4 ] [ 5 ] [ 6 ]. The work presented here focuses on processing sequencing data with the goal of delivering as many correct sequences as possible to the decoder. Similar work for genomic purpose has already been carried out, mainly in bioinformatics and metagenomics, but they mainly focuses on DNA clustering algorithms. Methods such as UCLUST [ 8 ], CD-HIT [ 9 ] and DNACLUST [ 10 ] use greedy algorithms for sequence similarity, while MeShClust [ 11 ] and MetaDEC [ 12 ] use machine learning approaches. Graph-based techniques, such as MMseqs2 [ 13 ], offer further strategies. However, a common limitation of these algorithms is their reliance on a fixed sequence identity threshold, which can be difficult to determine, particularly in DNA storage systems. Unlike genomic clustering, DNA data storage presents unique constraints due to shorter DNA strands (100-300 nucleotides). Although algorithms such as ALFATClust [ 14 ] and SEED [ 15 ] attempt to dynamically adjust clustering parameters, they remain optimized for metagenomic data. Recent advances in large-scale DNA clustering for storage applications include the work of Rashtchian et al. [ 16 ] and GradHC [ 17 ], both of which are capable of clustering billions of reads and producing highly reliable clustering results. The output of these programs are clusters containing a small number of reads, which are then used to build a consensus sequence. Our approach is distinct and aims to merge the clustering and consensus steps into a single process. To achieve this, we draw inspiration from bioinformatics techniques used in genome assembly, particularly the use of de Bruijn graphs [ 18 ] [ 19 ]. This method has already been explored by [ 20 ] and has proven to be particularly robust in handling DNA breaks, rearrangements, and indels. However, as the size of the encoded files increases, the construction of large de Bruijn graphs can lead to excessive computation times. Our method reduces algorithmic complexity by hierarchically dividing the initial set of reads into smaller subsets of predefined size. Each subset is expected to contain a portion of the original sequences along with their noisy copies. A de Bruijn graph is then constructed for each subset. Given the known length (N) of the target sequences and the nucleotide sequences that mark their start and end (the primers), the sequences to be reconstructed correspond to all paths of length N between two k-mers derived from the start and end primers. The advantages of this method are: (a) the hierarchical partitioning is extremely fast; (b) constructing smaller de Bruijn graphs I (c) parallelization on multicore architectures is straightforward since there are no dependencies between the de Bruijn graphs. 2 Method In the context of DNA storage, reference sequences consist of three concatenated parts, as illustrated in Figure 1 : a start primer, a payload, and an end primer. The payload carries the encoded data, along with indexing information for document reconstruction, optional checksums, and error correction data. To retrieve a document, a PCR is first performed using the specific primer pair associated with it, followed by a sequencing process. The resulting dataset may contain billions of reads. Here, we assume that each read includes a complete sequence, encompassing both primers. Download figure Open in new tab Figure 1: Structure of a reference sequence in the context of DNA storage To recover the initial set of sequences, a straightforward strategy would be to construct a de Bruijn graph using all these sequences. Since all sequences share the same start and end primers, one could then enumerate all paths in the graph that begin with a k-mer extracted from the start primer and end with a k-mer from the end primer. Each path represents a possible sequence. Paths that do not have the correct length are discarded. In practice, this solution works well and remains relatively fast for graphs of reasonable size. However, as the number of sequences increases, the computation time grows proportionally. Indeed, constructing a de Bruijn graph involves k-mer counting, which is a computationally intensive task. Also, identical k-mers in different sequences lead to an increase in the number of possible paths, drastically raising the complexity of the problem. To mitigate this complexity, our method first partitions the dataset (fastq file) into multiple subsets of approximately equal size. Each subset, and it is of the most importance, contains all instances of a given original sequence (see Section Partitioning). From there, as mentioned earlier, the problem consists of searching for paths in graphs that are now of reduced size. The algorithm represented figure 2 can be abstracted as follows: Download figure Open in new tab Figure 2: Overview of our method. Step 1 aligns the start primers and partition all reads. Step 2 builds de Bruijn graphs. Step 3 traverses the graphs to find consensus sequences of the expected length. Partition the dataset into subsets of approximate equal size For each subset - Create a de Bruijn graph - Generate consensus sequences from the graph To accelerate computation in a multithreaded environment, step 1 can benefit from intra-read parallelism, while step 2 allows for intra-subset parallelism. The following subections decribe in detail each of these steps. 2.1 Partitioning The objective of this step is to distribute the entire set of reads into subsets so that all copies of the same sequence are present in only one subset. The first task is to identify the start primer within the read, considering that the method is independent of the sequencing technology and that the primer is not necessarily located at the very beginning of the read. To achieve this, a local alignment is performed for each read against the start primer (and its reverse complement). However, some reads may be lost if no sufficiently good alignment is found. This is more likely to occur for reads with a high error rate. But removing them is actually beneficial, as they could otherwise compromise the accuracy of the consensus search. From the primer alignment position, a few bases immediately following the primer are used as a partitioning key. To account for the unknown data distribution, we use a dynamic key size. This approach ensures a more balanced load distribution and prevents some subsets from having an excessive number of reads while others remain nearly empty. We start with keys of 4 bases and then recursively split each subset into 64, assigning one subset to each possible value of the next 3 bases, thereby expanding the key. This process continues until each subset contains no more than a predefined number of reads. Ultimately, this step produces a set of FASTA files, each containing approximately the same number of reads. Each subset is written to disk and reads are trimmed to eliminate irrelevant parts: segments preceding the starting primer and segments after the expected length of the sequence are removed. This trimming can be done because the characteristics of the original reference are known, in particular the length of the oligonucleotides that encode a file. More generally, the partitioning strategy accept several start primers as large files are often indexed in several parts, each identified by a specific pair of primers. Partitioning is therefore initially based on these different start primers. This partitioning method is very fast and efficient. It systematically groups reads with similar start sequences. With a low error rate, the majority of copies of the same sequence will end up in the same subset. However, errors on the first bases of the read will inevitably lead to incorrect read assignment. But with current error rates and advances in sequencing technologies and base-calling algorithms, the number of misassigned reads remains fairly low, and does not disrupt processing in subsequent steps. 2.2 de Bruijn graph building De Bruijn graphs are primarily built based on k-mer counting. The advantage of considering only small subsets (fewer than 10,000 reads) is that the resulting graph will also be of limited size. The k-mer length can be optimally adjusted, enabling extremely simple and fast counting techniques that do not require complex data structures. In this case, a simple hash table is used and the values of the counters are capped to 255 to fit into 8-bit variables. Finally, we discard the k-mers that appear less than a user-defined threshold. This filtering eliminates sequencing errors, but also k-mers from reads assigned to an improper subset. 2.3 Generating consensus sequences The last stage consists of considering the k-mers previously found as a de Bruijn graph and traversing its paths to find consensus sequences. The exploration begins at the node representing the starting k-mer, which is actually the rightmost k-mer of the start primer. Then, the next possible four bases are considered recursively. If the traversal reaches the leftmost k-mer of the ending primer and the current payload has the exact expected size, the sequence is considered a consensus. If the traversal results in a sequence that exceeds the expected consensus length, the search on that particular path is terminated. Since this exploration has exponential complexity, due to cycles in the graph, several techniques are employed to reduce processing time. First, the graph is progressively pruned by removing nodes with no outgoing edges. Second, when a path is terminated, the current sequence is analyzed for repeated k-mers. If repetitions are detected, one edge is removed from the graph to break the loop. This optimization is lossy, meaning that it may prevent some valid consensus sequences from being found. However, eliminating loops significantly accelerates the search process. As in any data transmission context, some sequences are expected to be missing. To address this, efficient techniques such as error correction codes can be applied. 3 Evaluation setup 3.1 Hardware Hardware-wise, we use a server fitted with an Intel® Xeon® Silver 4215 CPU @ 2.5 processor (Skylake architecture, 16 hyper-threaded cores) and 256 of DDR4 @ 2.4 RAM. The server operates on Debian 10. 3.2 Software ConCluD ( C onsensus algorithm based on Clu stering for DNA storage D ecoding) is the name of the C++ software developed for the method described in the previous section. Unless state otherwise, ConCluD is executed with 32 CPU threads with the following default parameter values: k-mer size = 15 k-mer abundance threshold = 3 maximum number of reads in a subset = 5000 DBGPS has been developed by [ 20 ] and also uses a de Bruijn graph approach to reconstruct sequences after sequencing. A single de Bruijn graph is built from the whole read dataset. To reconstruct sequences, the oligonucleotides must be necessarily structured into three distinct fields: an index area, followed by the payload, and then a CRC field. The entire sequence is also framed by two primers. The C multithreaded version is used for comparison. GradeHC : has been developed by [ 17 ] and only generates clusters of reads. Each cluster can then be used to generate a single consensus sequence. 3.3 Datasets This section presents the different datasets we used for evaluating ConCluD. The following table summarizes the synthetic datasets. For each dataset, it gives its name, the size of the oligonucleotides, the sequencing technology, the number of primers, the number of oligonucleotides and the resulting reads with an average coverage of 30x. A detailed description of each of them, and the way they have been designed, can be found in annex. View this table: View inline View popup Download powerpoint The data are from high-resolution images. The payloads were constructed in compliance with the requirements for DNA storage (no large homopolymers, GC content around 50%, primers that do not hybridize in the middle of the payload, etc.). In addition, to ensure compatibility with the DBGPS dataset format, the same oligonucleotide structure was used for the IM1, IM10 and IM100 datasets, although ConCluD does not need any indexes to process the data. The ONT and Illumina sequencing datasets have been respectively generated with the PBSIM3 [ 22 ] and ART simulator tools [ 21 ]. The coverage has been approximately set to 30X. The three datasets IM1-xxx have been designed from small size JPEG images (497 KB and 1.3 MB) encoded with oligonucleotides of sizes 164nt and 272nt and synthetized with ONT or Illumina read simulators. Each image is encoded within a single pool of oligonucleotides sharing the same start and end primers, leading to 2 primers to perform a PCR. The 3 datasets IM10-xxx have been designed from larger images (5.8 MB and 13.5 MB) encoded with oligonucleotides of size 164 nt and 272 nt. Each image is encoded within 10 pools of oligonu-cleotides with 2 variants: for IM10-xxx20 each pool have a different pair of primers, leading to a total of 20 primers (10 × 2). For IM10-ont11, each pool have the same start primer, but a different end primer, reducing the number of primers to 11. The dataset IM100-272ont20 has been designed from a single very large JPEG image (134.7 MB) which is encoded into 100 different pools of oligonucleotides of size 272 nt. Each pool is distin-guished by a different pair of primers among 10 start and 10 end primers. The reason for splitting the binary file into small chunks of oligonucleotide pools is that it will make the decoding process simpler and faster in the downstream steps. Running error correcting codes on small pools, especially to retrieve lost oligonucleotides, is much less complex than doing it on large pools. In this experimentation it is also a way to investigate the best choice regarding pair of primers. We will see in the next section that sharing primers provides better performances. The datasets SF-xxx refer to the same image, but encoded with oligonucleotides of size 480nt, 840nt and 1740nt. The purpose of these three datasets is to test how ConCluD performs on large oligonucleotides. Since DBGPS primarily targets small oligonucleotides, these datasets won’t be used for comparison. To test and evaluate ConCluD on real data, the following real datasets have been selected: View this table: View inline View popup Download powerpoint The Park dataset was designed using Fountain codes to explore how low-quality reads can be leveraged to reduce sequencing coverage requirements [ 26 ]. The Goldman dataset [ 27 ] includes five files totaling 757 KB of data, in which text, images, and audio have been encoded. The Odeuropa dataset is provided by the OligoArchive European project [ 24 ], and the Dunhuang-dX datasets originate from [ 20 ] where ten digital pictures of Dunhuang murals have been encoded. Each Dunhuang datasets is based on the same initial sample but was subjected to accelerated aging for 0, 28, 56, or 70 days to induce molecular degradation. Finally, the HorrorMovie dataset is derived from Malcolm Le Grice’s Horror Film 1, which was encoded into nearly 2 million of nucleotides [ 23 ]. Illumina sequencing generates paired-end reads, which need to be merged to reconstruct complete oligonucleotides. This step is crucial for our strategy, as it requires both the start and end primers to be present on the same read. Merging was performed using the NGmerge tool [ 25 ]. However, in some datasets, primer or adapter sequences were missing due to trimming during the base calling process, which fortunately reorients all sequences identically. To maintain compatibility with our approach, artificial start and end primers have been added to restore the expected read structure. 3.4 Metrics To evaluate the accuracy of our algorithm, we consider the precision and recall metrics. Precision measures the amount of sequences returned by the program that are actually present in the reference. Achieving a good accuracy is important to ensure the reconstruction does not alter the original data. We compute the raw precision of the algorithm, independently of any checksum that might be embedded into the payload, and the final precision that allows to discard sequences after checksum verification. On the other hand, recall measures the proportion of the reference that was found by the program. A bad recall hurts the reconstruction of information as some data blocks may be missing. We compute the raw recall of the program, independently of any error correction codes or other device that might be embedded into the payload to recover missing sequences. Concerning time performance, we report the elapsed time of any execution. 4 Evaluation results 4.1 Precision and Recall Figure 3 reports the precision and recall of ConCluD and DBGPS software for each datasets. ConCluD parameters (k-mer size, k-mer threshold, subset size), have been set to optimized the results according to the dataset features. Overall, on small oligonucleotides (164nt and 272nt), ConCluD achieves excellent results in both precision and recall. On the other hand, DBGPS only produces satisfactory results for small datasets with limited oligonucleotide sizes. For the five last datasets, DBGPS does not output any consensus. More specifically, for the SF datasets, DBGPS crashes with a stack smashing error. For dataset IM10-272ont11, the execution time is excessively long (24 minutes) and the output consensus sequences are erroneous. And for dataset IM100-272ont20, the process was killed due to memory overflow (256 GB). Indeed, in DBGPS, the size of the de Bruijn graph is directly proportional to the number of reads (and thus to the number of k-mers). The complexity of traversing the graph and reconstructing paths increases accordingly. Removing loops in the graph to reduce complexity has an immediate impact: many paths supporting reference sequences are not explored. ConCluD integrates the same strategy but applies it to a multitude of small graphs, where only a few loops are removed, leading to marginal impact on the results. Download figure Open in new tab Figure 3: Precision and Recall on the synthetic datasets. Raw precision measures the amount of sequences returned by the software without CRC. Overall precision, actually 100% on each datasets, take care of CRC As previously mentioned, the oligonucleotides were structured to be compatible with the DBGPS format. ConCluD relies on the nucleotide sequence directly following the starting primer, which, in this case, corresponds to the added index and could introduce bias into our strategy. Therefore, we tested our partitioning technique by using the partitioning position after the index, which corresponds to the beginning of the payload. The results remain identical. No bias is introduced. The last three datasets (SF-xxx) test the robustness of ConCluD against long oligonucleotides. A slight decrease in recall is observed when the length of the oligonucleotides increases. On the raw precision side, however, a very clear decrease is observed, which is consistently compensated by a validity check. This evaluation shows that achieving 100% recall is never possible, either because on large oligo pools the synthesis did not produce all the expected oligonucleotides or because the random sequencing process failed to capture all the information. This is why error-correcting codes are added to recover the information. Although this is not the focus of this work, these results may help quantifying sequencing losses and adjusting the parameters of error-correcting codes. 4.2 Execution time Figure 4 shows the execution times for ConCluD and DBGPS. For the datasets IM10-272ont11 and IM100-272ont20, DBGPS execution times are not reported as no results are produced. For ConCluD, two things can be observed: Download figure Open in new tab Figure 4: Execution time for ConCluD and DBGPS. The execution time increases almost linearly with the data size. The top of the figure shows the relative execution time using the smallest dataset execution time as reference. The difference between IM10-272ont11 (with a common start primer for all images) and IM10-272ont22 (with different start and end primers for each image) favors the dataset using only 11 primers. There is a significant decrease, which can be explained by the optimization of the graph traversal. Instead of performing 10 searches with 10 different start primers, only one search is initiated. What’s more, instead of performing 10 alignments for each primer, we do just one, so partitioning also involves 10x fewer calculations. Moreover, it’s better to have 1 start primer and 10 end primers, than 10 start primers and 1 end primer. Behavior is not symmetrical. GradHC, that only performs clustering for DNA storage, has been run on dataset IM1-164ont2 (1 image - 31,294 oligonucleotides). The measured execution time was 2580 seconds. This time is actually far too long for this method to be practically used. ConCluD has two major steps: (1) partitioning the dataset; (2) building consensus sequences from a de Bruijn graph. The table below indicates the execution time (in second) for each step: View this table: View inline View popup Download powerpoint We immediately notice that the partitioning step is dominant. The construction of de Bruijn graphs (mainly k-mer counting) is fast, as well as the path search within these different graphs. 4.3 Parameter impacts There are three main parameters that govern the operation of ConCluD: the size of the subsets in the partitioning phase, the size of the k-mers and their solidity in the consensus phase. These parameters can influence the quality of the results (precision and recall) as well as the execution time. Experimental tests show that a subset size of 5,000 reads, a k-mer size of 15 nt, and a low-abundance k-mer threshold of 3 provide good results and can be used as default parameters. In the following sections, the impact of these parameters is evaluated separately, with the others set to their default values. Figure 5 evaluates the impact of subset size using three different values: 5K, 25K, and 50K reads. The results clearly show that larger subsets significantly reduce recall. This supports our strategy of managing multiple small de Bruijn graphs rather than a single large one. The execution time is less affected; partitioning into larger batches is slightly faster. However, the processing time for consensus finding increases significantly. A maximum subset size of 5K reads therefore appears to be an optimal choice. Download figure Open in new tab Figure 5: Impact of the subset size for precision, recall and execution time The size of the k-mers ( figure 6 ) have a strong impact. As it can be seen the minimum size is 13nt. Above this value, both precision and recall are good, and the computation time remains stable. A default value of 15 is therefore entirely appropriate. Download figure Open in new tab Figure 6: Impact of the k-mer size for precision, recall and consensus execution time If the value of the low-abundance k-mer threshold ( Figure 7 ) used to build the de Bruijn graph is too low (i.e. set to 1) too many false positive sequences are generated, leading to very poor precision and excessive computation time due to very noisy de Bruijn graph. On the other hand, as this value increases, the recall decreases accordingly. A good balance to optimize the recall is achieved with a trade-off value of 3 for an expected coverage of 30X. Download figure Open in new tab Figure 7: Impact of the k-mer abundance threshold for precision, recall and consensus execution time 4.4 Evaluation on real datasets Finally, ConCluD was run and evaluated on real datasets, which were subsampled when required to ensure a maximum coverage of 30X. The results are presented in the following table. View this table: View inline View popup Download powerpoint For the Dunhuang-dX datasets, we report only the raw precision, as the verification protocol is unavailable. However, when this information is accessible, incorporating it typically increases the precision to 100%. Recall remains high, except for degraded sequences in the Dunhuang datasets. In such cases, our strategy becomes ineffective, as the molecular damage fragments the sequences, preventing ConCluD from performing consistent partitioning. The experiments conducted by the authors of DBGPS represent extreme conditions, and it is reasonable to assume that a well-designed archival system would avoid such scenarios. 5 Conclusion We presented a method based on de Bruijn graphs to process sequencing reads derived from oligonucleotides and reconstruct the DNA sequences encoding a digital document. The method has been implemented in the ConCluD software. The originality and efficiency of our approach lie in a partitioning strategy that groups each sequence and all of its copies into the same subset. Each subset is then processed independently. This offers two key advantages: first, the resulting graphs are smaller, which simplifies their construction and handling; second, this execution model allows for highly scalable processing. Our strategy is also compatible with longer oligonucleotide lengths, anticipating future technological advances that could enable the storage and retrieval of significantly larger molecules than what is currently possible. Our method does not require any prior knowledge of the payload’s structure. However, it does assume that each read contains a complete oligonucleotide, including its start and end primers. These primers are essential to initiate and terminate path-finding in the de Bruijn graphs. Unfortunately, some Illumina paired-end sequencing protocols do not retain this information. The base calling step typically removes primers, as they are biologically irrelevant in that context. Nevertheless, the oligonucleotide can be reconstructed by merging read pairs and appending artificial primers. As such, the ConClud software is compatible with both long- and short-read sequencing technologies. DNA-based data storage is a promising solution for long-term, high-density information archiving. In this context, recovering high-quality sequences after the sequencing step is of utmost importance, as it directly influences the design of error-correcting codes tailored for DNA storage. When sequence reconstruction is highly accurate, the level of redundancy required for error correction can be significantly reduced, leading to substantial cost savings in the synthesis process. Moreover, being able to efficiently process large volumes of sequencing data is essential to ensure a smooth and scalable end-to-end data storage pipeline. Our work contributes to addressing both of these challenges, offering a reliable and efficient solution for DNA data storage systems. Acknowledgments The work described in this paper is funded by the French ANR program ANR-22-PEXM-0005 MoleculArxiv and by the Labex CominLabs (dnarXiv project). Annex 1: ConCluD Software Two versions of ConCluD exist: one designed for conventional multicore servers and one adapted to Processing-in-Memory architectures. Here, only the multicore version is considered. The source code is available here: https://gitlab.inria.fr/pim/org.pim.dnarxiv/-/tree/main Annex 2: Generation of the IM1, IM10 and IM100 synthetic datasets The datasets used for performance evaluation have been constructed from datasets archived in the following Zenodo repository: https://zenodo.org/records/15387164 6 datasets that represents DNA encoded JPEG images are available: IM1-140-2 : 1 pool of 31,294 oligonucleotides of size 140 nt. All oligonucleotides share the same pair of primers (2 primers). IM1-248-2 : 1 pool of 35,090 oligonucleotides of size 248 nt. All oligonucleotides share the same pair of primers (2 primers). IM10-140-20 : 10 pools of 36,583 oligonucleotides of size 140 nt. Each pool have a different pair of primers (20 primers). IM10-248-20 : 10 pools of 35,069 oligonucleotides of size 248 nt. Each pool have a different pair of primers (20 primers). IM10-248-11 : 10 pools of 35,069 oligonucleotides of size 248 nt. Each pool have the same start primer but a different end primer (11 primers). IM100-248-20 : 100 pools of 35,114 oligonucleotides of size 248 nt. Each pool have a different pair of primers among 10 start primers and 10 end primers (20 primers). The structure of an oligonucleotide is as follows: Download figure Open in new tab The length of the primers are 20 nt. The length of the payload is 98 nt for oligos of size 140 nt and 206 nt for oligos of size 148 nt. A 1-bit buffer is inserted between the payload and the primers due to biotechnological constraints. The payload contains encoded indexing and CRC information. The scripts for generating the oligonucleotide datasets are available here: https://gitlab.inria.fr/pim/org.pim.dnarxiv/-/tree/main/paper_scripts/datasets/IM-1-10-100 The generation of N pools of oligonucleotides from a binary file is performed as follows: Split the binary file into N binary strings. N is determined by the target length of the oligonucleotides and the number of oligonucleotides by pool. Apply a SHA256 hash to each binary strings with known key in order to obtain (pseudo) random binary strings of identical size, which is added to the binary string with an XOR operation. This operation ensures a good GC percentage balance and avoids large homopolymers. Cut the binary strings into P fragments, relatively to the length of the oligonucleotides Append a 16-bit index and 20-bit checksum to each fragment. Convert fragments to DNA sequences with an encoding algorithm at a density of 1.66 bit/nt, ensuring a good GC balance and no homopolymeres. Compute primers that do not hybridize with the fragments of a given pool. Add primers and buffer to each fragment. In order to be compliant with the DBGPS software, oligonucleotides from these datasets are modified. A 4 bytes index (= 16 bases) is inserted between the start primer and the payload, and a 2 bytes CRC (= 8 bases) is inserted between the payload and the end primer. The index is chosen to not have any homopolymers larger than 4 bases. The CRC is computed with the code of the public Python implementation of DBGPS. The helper scripts to convert the reference prints in the console the length and index start/end parameters to use to run ConCluD and DBGPS. Long reads are generated with the PBSIM3 simulator, and short reads by the ART simulator. As a result we get the following synthetic datasets: IM1-164ont2 made from IM1-140-2 and PBSIM3 IM1-164ill2 made from IM1-140-2 and ART IM1-272ont2 made from IM1-248-2 and PBSIM3 IM10-164ont20 made from IM10-140-20 and PBSIM3 IM10-272ont20 made from IM10-248-20 and PBSIM3 IM10-272ont11 made from IM10-248-11 and PBSIM3 IM100-272ont20 made from IM100-248-20 and PBSIM3 Annex 3: Experiment Information All the conditions and the implementation of the various experiments carried out within the frame-work of this article are available here: https://gitlab.inria.fr/pim/org.pim.dnarxiv/-/tree/main/paper_scripts Funder Information Declared French ANR program MoleculArxiv , ANR-22-PEXM-0005 Labex CominLabs , dnarxiv Footnotes florestan.de-moor{at}irisa.fr olivier.boulle{at}irisa.fr lavenier{at}irisa.fr https://zenodo.org/records/15387164 References [1]. ↵ G. M. Church et al. , Next-Generation Digital Information Storage , DNA. Science , 337 , 2012 . [2]. ↵ L. Ceze , J. Nivala , K. Strauss , Molecular digital data storage using DNA , Nat. Rev. Genet . 20 , 456 – 466 , 2019 . OpenUrl CrossRef PubMed [3]. ↵ L.C. Meiser et al. , Reading and writing digital data in DNA , Nat. Protoc . 15 , 86 – 101 , 2020 . OpenUrl CrossRef PubMed [4]. ↵ W.H. Press , J.A. Hawkins , S.K. Jones , J.M. Schaub , I.J. Finkelstein , HEDGES error-correcting code for DNA storage corrects indels and allows sequence constraints , Proc. Natl. Acad. Sci. USA 117 , 18489 – 18496 , 2020 . OpenUrl Abstract / FREE Full Text [5]. ↵ Y. Erlich , D. Zielinski , DNA Fountain enables a robust and efficient storage architecture , Sci . 355 , 950 – 954 , 2017 . OpenUrl Abstract / FREE Full Text [6]. ↵ E. Marinelli , Y. Yan , L. Tattini , V. Magnone , P. Barbry , R. Appuswamy , CMOSS: A Reliable, Motif-based Columnar Molecular Storage System . In Proceedings of the 17th ACM International Systems and Storage Conference (SYSTOR ‘24) . Association for Computing Machinery , New York, NY, USA , 2024 [7]. ↵ H.M. Kiah , G.J. Puleo , O. Milenkovic , Codes for DNA Sequence Profiles , IEEE Trans. Inform. Theory 62 , 3125 – 3146 , 2016 . OpenUrl CrossRef [8]. ↵ R.C. Edgar , Search and clustering orders of magnitude faster than BLAST , Bioinformatics , 26 : 2460 – 1 , 2010 . OpenUrl CrossRef PubMed Web of Science [9]. ↵ L. Fu , B. Niu , Z. Zhu Z et al. , CD-HIT: accelerated for clustering the next-generation sequencing data , Bioinformatics , 28 : 3150 – 2 , 2012 OpenUrl CrossRef PubMed Web of Science [10]. ↵ M. Ghodsi , B. Liu , M. Pop , Dnaclust: accurate and efficient clustering of phylogenetic marker genes , BMC Bioinformatics , 12 : 271 , 2011 . OpenUrl CrossRef PubMed [11]. ↵ B.T. James , B.B. Luczak , H.Z. Girgis , MeShClust: an intelligent tool for clustering DNA sequences , Nucleic Acids Res , 46 : e83 , 2018 . OpenUrl CrossRef PubMed [12]. ↵ H.Q. Bao , L.V. Vinh , T. Van Hoai , A deep embedded clustering algorithm for the binning of metagenomic sequences , IEEE Access , 10 : 54348 – 57 , 2022 . OpenUrl CrossRef [13]. ↵ M. Steinegger , J. Söding , Clustering huge protein sequence sets in linear time , Nat Commun , 9 : 2542 , 2018 . OpenUrl CrossRef PubMed [14]. ↵ J.K.H Chiu , R.T. Ong , Clustering biological sequences with dynamic sequence similarity threshold , BMC Bioinformatics , 23 : 108 , 2022 . OpenUrl CrossRef PubMed [15]. ↵ E. Bao , T. Jiang , I. Kaloshian et al. , SEED: efficient clustering of next-generation sequences , Bioinformatics , 27 : 2502 – 9 , 2011 OpenUrl CrossRef PubMed Web of Science [16]. ↵ C. Rashtchian , K. Makarychev , M. Rácz , S.D. Ang , D. Jevdjic , S. Yekhanin , L. Ceze , K. Strauss , Clustering billions of reads for DNA Data storage , Advances in Neural Information Processing Systems , 2017 . [17]. ↵ D. B. Shabat , A. Hadad , A. Boruchovsky , E. Yaakobi , GradHC: highly reliable gradual hash-based clustering for DNA storage systems , Bioinformatics , Volume 40 , Issue 5 , May 2024 . [18]. ↵ D.R. Zerbino , E. Birney , Velvet: algorithms for de novo short read assembly using de Bruijn graphs . Genome Res . 18 , 821 – 829 , 2008 . OpenUrl Abstract / FREE Full Text [19]. ↵ P.A. Pevzner , H. Tang , M.S. Waterman , An Eulerian path approach to DNA fragment assembly , Proc. Natl. Acad. Sci. USA 98 , 9748 – 9753 , 2001 . OpenUrl Abstract / FREE Full Text [20]. ↵ L. Song , F. Geng , ZY Gong et al. , Robust data storage in DNA by de Bruijn graph-based de novo strand assembly , Nature Communications 13 , 5361 , 2022 . OpenUrl CrossRef PubMed [21]. ↵ W. Huang , L. Li , J.R. Myers , G.T. Marth , ART: a next-generation sequencing read simulator . Bioinformatics , 15; 28 ( 4 ), 2012 . [22]. ↵ Y. Ono , M. Hamada , K. Asai , PBSIM3: a simulator for all types of PacBio and ONT long reads , NAR Genomics and Bioinformatics , Volume 4 , Issue 4 , 2022 . [23]. ↵ L. Curham , L. Ihlein , R. Appuswamy , Re-enactment, users manuals and DNA storage: methods for media art , ISEA 2024, 29th International Symposium on Electronic Art , 21-29 June 2024 , Meanjin, Australia [24]. ↵ R. Appuswamy , K. Lebrigand , P. Barbry , M. Antonini , O. Madderson , et al. , OligoArchive: Using DNA in the DBMS storage hierarchy . Biennal Conference on Innovative Data Systems Research (CIDR 2019) , Jan 2019 , Asilomar, CA , United States . ffhal-03661386f [25]. ↵ JM. Gaspar , NGmerge: merging paired-end reads via novel empirically-derived models of sequencing errors . BMC Bioinformatics 19 , 536 , 2018 . [26]. ↵ J. Park et al. , Sequence analysis and decoding with extra low-quality reads for DNA data storage , https://github.com/PParkJy/SAD-DNAstorage [27]. ↵ N. Goldman , P. Bertone , S. Chenmet al ., Towards practical, high-capacity, low-maintenance information storage in synthesized DNA . Nature 494 , 77 – 80 , 2013 . OpenUrl CrossRef PubMed Web of Science View the discussion thread. Back to top Previous Next Posted May 23, 2025. Download PDF Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following De Bruijn Graph Partitioning for Scalable and Accurate DNA Storage Processing Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share De Bruijn Graph Partitioning for Scalable and Accurate DNA Storage Processing Florestan De Moor , Olivier Boullé , Dominique Lavenier bioRxiv 2025.05.19.654814; doi: https://doi.org/10.1101/2025.05.19.654814 Share This Article: Copy Citation Tools De Bruijn Graph Partitioning for Scalable and Accurate DNA Storage Processing Florestan De Moor , Olivier Boullé , Dominique Lavenier bioRxiv 2025.05.19.654814; doi: https://doi.org/10.1101/2025.05.19.654814 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Bioinformatics Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41913) Biophysics (21436) Cancer Biology (18578) Cell Biology (25482) Clinical Trials (138) Developmental Biology (13372) Ecology (19889) Epidemiology (2067) Evolutionary Biology (24302) Genetics (15599) Genomics (22483) Immunology (17728) Microbiology (40365) Molecular Biology (17163) Neuroscience (88540) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4821) Physiology (7637) Plant Biology (15130) Scientific Communication and Education (2045) Synthetic Biology (4290) Systems Biology (9818) Zoology (2269)