Sama Algorithm Introduces Correctness Guarantees for De Bruijn Graph Genome Assembly

A computational framework and genome assembly tool named Sequence Assembly avoiding MisAssemblies introduces a mathematical model to calculate the probability of misassembly at each position of a de Bruijn graph based assembly while explicitly accounting for missing data. Developed by Leena Salmela and published in Algorithms for Molecular Biology, the tool establishes upper bounds on structural assembly errors and generates contigs accompanied by position-specific correctness estimates.

Genome assembly in computational biology addresses the problem of reconstructing a target genome from sequencing reads, historically formulated through shortest common superstrings, Hamiltonian paths in string graphs, or Eulerian tours in de Bruijn graphs. Because practical solutions to these formulations are almost never unique, contemporary assemblers rely on heuristic procedures that remain difficult to analyze and frequently introduce structural misassemblies into output sequences.

Current assembly algorithms generally output unitigs, which correspond to nonbranching paths in string or de Bruijn graphs that are guaranteed to represent substrings of the genome only in the complete absence of missing graph branches. More recent theoretical formulations such as omnitigs guarantee valid substrings across all genomes compatible with the graph, but they likewise assume complete data coverage, causing both paradigms to risk reporting incorrect sequences when missing data occurs.

Previous probabilistic assembly methods have explored maximum likelihood reconstructions, simulated annealing schemes, or Markov chain Monte Carlo approaches, yet these frameworks either did not calculate positional misassembly probabilities, failed to address missing data, or proved too computationally demanding for organisms larger than bacteriophage PhiX174. Other investigations established that polynomial-time unique reconstructions require all repeats to be spanned by reads, an assumption that cannot accommodate missing data or supply localized confidence estimates.

To overcome these limitations, the newly presented method constructs a de Bruijn graph from read sets by cataloging k-mers as nodes and connecting nodes that overlap by k minus 1 characters with edges representing k plus 1-mers. For every node and edge, the algorithm tallies read abundance to evaluate whether an edge representing a one-character extension corresponds to a unique continuation in the underlying genome or poses a misassembly risk.

For any sequence S occurring a times in the read set and possessing an extension Sc, the probability of Sc not being a unique extension is evaluated based on whether S represents a repeat with alpha occurrences. In the worst-case configuration where alpha minus 1 genomic extensions are identical and one differs, the observed abundance of Sc follows a binomial distribution with parameters a and (alpha minus 1) divided by alpha, which requires edge abundances to exceed ((alpha minus 1) divided by alpha) times a.

By applying a right-tail upper bound established by Ferrante, the model bounds misassembly probabilities using odds ratios and Kullback-Leibler divergences between Bernoulli random variables. Combining these bounds with repeat multiplicity proportions obtained via the Detox tool, which estimates repeat distributions for alpha values between zero and five from k-mer abundance histograms, the algorithm calculates total probabilities of non-unique extensions across all repeat classes.

In practical implementation, the pipeline constructs de Bruijn graphs using BCALM 2 and establishes abundance thresholds across all k-mer frequencies for a chosen misassembly threshold epsilon, avoiding computationally expensive per-edge evaluations. The contig generation algorithm prunes nodes with abundances below five to remove sequencing errors and traverses maximal paths exclusively along edges that satisfy abundance thresholds bidirectionally, resolving potential conflicts at repeat boundaries and error junctions where unidirectional extensions occur.

The assembly tool was evaluated on real PacBio HiFi reads of Escherichia coli downsampled to 20x, 40x, and 80x coverage, real Illumina reads of Staphylococcus aureus downsampled to 80x coverage, and simulated human chromosome 21 read sets at 40x and 60x coverage generated using Art. Benchmark tests compared the tool against unitigs from BCALM 2 and contigs from SPAdes without repeat resolution across k-mer lengths of 31 and 63 alongside misassembly probability thresholds ranging from 0.4 down to 0.0001.

Evaluation with QUAST demonstrated that higher coverage datasets allowed the tool to achieve contiguity metrics, such as NGA50 and genome fraction, comparable to heuristic tip and bulge removal in SPAdes while completely avoiding misassemblies. Across all benchmark datasets, setting the misassembly threshold to 0.01 or lower eliminated misassemblies entirely, whereas misassemblies only emerged on the more complex human chromosome 21 datasets under relaxed probability thresholds.

Runtime and memory evaluations on an eight-core machine confirmed that the pipeline remains competitive with existing tools even when accounting for the execution of BCALM 2 and Detox. The author noted that while the current model is limited to single haploid genomes, its position-level probability bounds provide a mathematical basis for evaluating structural variant confidence, weighting de Bruijn graph sequence alignments, and extending future assembly guarantees to diploid genomes and metagenomics.

In de Bruijn graph representations of sequencing data, an assembly corresponds to a collection of directed paths connecting overlapping k-mers. Each path sequence begins with the full k-character string of the initial node and appends the final nucleotide of each subsequent node in the path. Within this graph structure, establishing structural correctness requires determining whether an edge representing a one-base extension Sc of a k-mer S reflects an unambiguous continuation in the target genome. If a sequence S occurs only once in the genome, the probability that an observed extension is non-unique equals zero. Even when a k-mer represents a repeat sequence, an extension remains uniquely valid if every genomic copy of S is followed by the identical nucleotide, establishing equal multiplicity between the k-mer and its extension. The core computational challenge, formulated as Problem 1, centers on calculating the probability that an extension Sc attains an observed read abundance of at least ell when S occurs as a genomic repeat possessing two or more distinct extensions.

To establish bounds on this probability, the mathematical analysis models an alpha-repeat where S appears alpha times across the genome, yielding alpha genomic extensions. The most difficult scenario to resolve arises when alpha minus one extensions share an identical nucleotide while exactly one extension diverges. Under this worst-case configuration, the observed frequency of the dominant extension follows a binomial distribution over the total k-mer abundance a, requiring the extension abundance threshold to exceed the fraction of alpha minus one divided by alpha multiplied by a. Applying Ferrante’s right-tail inequality bounds the cumulative tail probability using the odds ratio and the Kullback-Leibler divergence between Bernoulli random variables. Because any of the alpha extensions could manifest this asymmetric distribution, multiplying by alpha bounds the likelihood that any extension exceeds the threshold. The framework integrates these values across repeat classes, weighting each term by the genomic proportion of alpha-repeats estimated from k-mer frequency histograms.

Rather than recomputing these probability integrals for every edge in the graph, the implementation predetermines minimum abundance thresholds for each k-mer abundance level based on the chosen error tolerance epsilon. A pivotal property of these derived thresholds is that they consistently exceed half of the originating k-mer abundance. As a result, at most one outgoing edge from any given node can satisfy the abundance requirement, which inherently prevents branching ambiguity and eliminates the need to resolve complex graph forks. However, directional discrepancies often occur at repeat borders and sequencing error boundaries. At the boundary of a repeat, a non-repeat k-mer can uniquely extend toward a repeat node, whereas the reverse extension cannot be verified as unique. Similarly, an erroneous k-mer absent from the genome can uniquely extend into a genuine k-mer, but the reverse transition remains ambiguous. By enforcing that edges must represent unique extensions in both forward and reverse directions, the traversal algorithm suppresses these false connections while generating contigs through bidirectional traversals from nodes with abundance values of five or greater.

The empirical validation examined specific dataset characteristics and parameter interactions. Choosing odd values of k, specifically 31 and 63, fulfilled algorithmic requirements while enabling compact data packing into 64-bit and 128-bit machine words. Across these settings, increasing k steadily enhanced assembly continuity. Threshold curves plotted across experimental read sets revealed distinct behavior between coverage levels. For k-mers with low abundance, thresholds clustered near half the k-mer frequency because unique sequences carry zero theoretical misassembly risk. For high-abundance k-mers, required thresholds approached full k-mer abundance, reflecting the heavier evidence burden needed to confirm unique extensions among repetitive sequences. These curves progressed unevenly because the proportion of repeat copies shifted across abundance levels. Furthermore, under highly stringent misassembly probabilities, certain k-mer abundances lacked any valid threshold capable of guaranteeing the requested error limit. The analysis also highlighted coverage nuances in the Staphylococcus aureus benchmark: despite an 80x nucleotide coverage, short 101-base-pair read lengths reduced effective k-mer coverage to 56x for 31-mers and 31x for 63-mers.

In the comparative evaluations, SAMA delivered higher contiguity than BCALM 2 across appropriate misassembly thresholds. When contrasted with SPAdes, SPAdes achieved somewhat higher contiguity on datasets with reduced k-mer coverage, such as the 20x Escherichia coli and Staphylococcus aureus sets, whereas SAMA achieved superior contiguity on the more structurally complex human chromosome 21 benchmarks. The author noted that resource comparisons against SPAdes were influenced by unavoidable repeat resolution routines in SPAdes. The implementation, distributed under the GNU AGPL license on GitHub, is coded in C++ for Linux systems, utilizing CMake, GCC, Google sparsehash, and GNU Multiple Precision libraries alongside Boost for approximate probabilistic inference.

Subject of Research: Structural correctness guarantees and misassembly probability estimation in de Bruijn graph genome assembly

Article Title: Sama: a contig assembler with correctness guarantee

Article References: Salmela, L. (2025). Sama: a contig assembler with correctness guarantee. Algorithms for Molecular Biology, 20(1), Article 9. https://doi.org/10.1186/s13015-025-00280-y

Image Credits: AI Generated

DOI: 10.1186/s13015-025-00280-y

Keywords: genome assembly, de Bruijn graph, contig generation, misassembly bounds, bioinformatics algorithms, computational biology, repeat resolution, sequencing coverage, Sama, contig, assembler, correctness