De Bruijn Graph Assembly, in depth: k-mers, solid filtering, unitigs, tips, bubbles, repeats and choosing k

From reads to a graph of k-mers

Take k = 5 and the read ATCGGATT. Its k-mers are ATCGG, TCGGA, CGGAT, GGATT. Each k-mer becomes a directed edge from its first four bases to its last four, so ATCGG is an edge from ATCG to TCGG. Consecutive k-mers in a read share four bases, so the read becomes a path through the graph. Two reads that overlap by at least k bases share k-mers, and so share edges, without anyone having compared them.

Building the graph needs a single pass that counts every k-mer in every read. Construction is linear in the total number of bases, and the graph has at most one node per distinct (k-1)-mer, whatever the read count. Deeper coverage raises the counts but does not add nodes, which is why de Bruijn assemblers scale to datasets that overlap-based assemblers could not handle.

Reverse complements. DNA is double-stranded, and a read is equally likely to come from either strand. The k-mer CGGAT on one strand is ATCCG on the other: reverse the string and swap A with T and C with G. Assemblers therefore store a canonical k-mer, the lexicographically smaller of a k-mer and its reverse complement, and add both orientations as edges when building the graph. Odd k is standard because an odd-length string can never equal its own reverse complement. The middle base would have to be its own complement, and no base is. That removes a class of self-loops that would otherwise confuse the traversal.

COMP = str.maketrans("ACGT", "TGCA")

def revcomp(s):
    return s.translate(COMP)[::-1]

def canonical(kmer):
    rc = revcomp(kmer)
    return kmer if kmer <= rc else rc

def count_kmers(reads, k, both_strands=True):
    counts = Counter()
    for r in reads:
        for i in range(len(r) - k + 1):
            km = r[i:i + k]
            if "N" in km:                 # ambiguous base: skip the k-mer
                continue
            counts[canonical(km) if both_strands else km] += 1
    return counts

Worked example: assembling 18 bases

The genome is ATCGGATTACCGGATGCC, 18 bases, containing the 5-mer CGGAT twice. To keep the trace readable this example uses one strand only. The reads are every 8-base window starting at an even position, sampled three times, plus one extra read and one read with an error, the C at position 9 (counting from 0) read as G: 20 reads in all. Counting 5-mers gives this histogram of how many distinct k-mers occur how often:

The three k-mers seen once are ATTAG, TTAGC and TAGCG, exactly the k-mers that overlap the error. This is the central fact of de Bruijn assembly: one substitution corrupts up to k k-mers, but those k-mers are rare, while true k-mers recur at roughly the coverage. The k-mer seen 13 times is CGGAT, the repeat. Then build the graph and compact every non-branching path into a unitig:

def build_graph(kmers, both_strands=False):
    out, inn = defaultdict(list), defaultdict(list)
    for km in kmers:
        for e in ((km, revcomp(km)) if both_strands else (km,)):
            u, v = e[:-1], e[1:]          # (k-1)-mer prefix -> suffix
            out[u].append(v)
            inn[v].append(u)
    return out, inn

def unitigs(out, inn):
    simple = lambda n: len(out[n]) == 1 and len(inn[n]) == 1
    result = []
    for u in list(set(out) | set(inn)):
        if simple(u):
            continue                       # unitigs start at branching or end nodes
        for v in out[u]:
            path = u + v[-1]
            while simple(v):               # extend through 1-in-1-out nodes
                v = out[v][0]
                path += v[-1]
            result.append(path)
    return result                          # a graph that is one pure cycle needs a
                                           # separate pass; omitted here

counts = count_kmers(reads, 5, both_strands=False)   # one-strand toy data
solid = [km for km, c in counts.items() if c >= threshold]
print(sorted(unitigs(*build_graph(solid))))

For real reads, pass both_strands=True to both functions; every contig then also appears as its reverse complement. With threshold 1 the run prints six unitigs: ATCGGA, ATTACCGGA, ATTAGCG, CGGAT, GGATGCC and GGATTA. The error has split the loop at node ATTA and hung a dead end, ATTAGCG, off it. With threshold 3 the three error k-mers vanish and four unitigs remain: ATCGGA, CGGAT, GGATTACCGGA and GGATGCC. The graph in the figure has exactly one Eulerian walk: start, the repeat, the loop, the repeat again, the exit. It spells the genome. Notice how the repeat appears as a single edge with double coverage rather than as two copies. Without coverage information, or a path that must use every edge, nothing in the graph says how many times to traverse it.

Cleaning the graph: errors, tips and bubbles

Real data has error rates around 0.1 to 1 percent on short-read platforms, and errors produce three recognisable shapes in the graph. Assemblers remove them in roughly this order.

Solid k-mer filtering. Plot the k-mer count histogram before anything else. A well-behaved library shows a large spike at count 1 and 2 (errors), a trough, and a broad peak at the k-mer coverage (genomic k-mers), with smaller peaks at multiples of it for repeats and, in diploid samples, a peak at half coverage for heterozygous k-mers. Set the threshold in the trough. Counting tools such as Jellyfish and KMC exist mainly to produce this histogram quickly.

Tips. An error near the end of a read leaves a short dead-end branch, as in the example. Remove a branch that ends with no successor if it is shorter than about 2k edges and its coverage is well below the branch it leaves.

Bubbles. An error in the middle of a read creates an alternative path of about k edges that leaves the true path and rejoins it. Velvet's tour bus algorithm and later variants detect two paths between the same pair of nodes with similar length, and merge the lower-coverage one into the higher. Heterozygous variants produce bubbles too, with balanced coverage, so popping them collapses the two haplotypes into one sequence. That is usually what a reference assembly wants, and it is exactly wrong for variant calling.

Chimeric and low-coverage links. Edges whose coverage is far below both neighbours are cut. This is the most dangerous step, because genuinely low-coverage regions, such as GC-rich sequence that amplifies poorly, look the same.

def k_mer_coverage(read_cov, read_len, k):
    # every read of length L contributes L - k + 1 k-mers
    return read_cov * (read_len - k + 1) / read_len

def p_kmer_error_free(err_rate, k):
    return (1 - err_rate) ** k

# 30x of 150 bp reads:  k=31 -> 24.0x k-mer coverage, 73% of k-mers error-free at 1% error
#                       k=55 -> 19.2x,                  57% error-free

Repeats, and why the output is contigs

A repeat shorter than k-1 bases is invisible: no (k-1)-mer occurs twice, so no node branches. A repeat longer than k collapses into one unitig with several entrances and exits, as CGGAT did. The graph then admits many Eulerian walks, and only one of them is the genome. Real genomes contain transposons thousands of bases long and segmental duplications far longer, so the Eulerian path is ambiguous almost everywhere at the chromosome scale.

This is why practical assemblers do not output an Eulerian path. They output contigs: unitigs, or paths through regions that the evidence resolves unambiguously, and stop at every unresolved branch. Three kinds of evidence extend them. Paired-end reads, two reads from the ends of a fragment of known approximate length, show which exit of a repeat follows which entrance. Coverage tells how many copies a collapsed repeat stands for. Long reads that span the repeat settle it directly, which is why hybrid and long-read pipelines produce far more contiguous assemblies.

Choosing k

k sets a trade-off with no free lunch. A larger k resolves more repeats, since any repeat shorter than k-1 disappears, but it lowers k-mer coverage and raises the chance that a k-mer contains an error, so low-coverage regions fragment into gaps. A smaller k connects through low coverage and errors but collapses more repeats. The formulas above make it concrete: on 30x of 150 bp reads with 1 percent errors, moving k from 31 to 55 cuts k-mer coverage from 24x to 19x and the error-free fraction of k-mers from about 73 percent to 57 percent.

Modern assemblers avoid choosing a single value. SPAdes and MEGAHIT iterate over several k values, building the graph for a small k, using its contigs as extra input for the next larger k, and keeping connectivity from the small graph and repeat resolution from the large one.

Memory at genome scale

The human genome has about 3.1 billion bases and roughly as many distinct genomic k-mers for k around 31. At 30x coverage with 1 percent error, the number of distinct erroneous k-mers can be several times that, since almost every error produces new ones. A k-mer up to 32 bases packs into one 64-bit word at 2 bits per base, but a hash table with counts and load-factor slack still costs well over 10 bytes per entry, so naive counting of a human sample needs hundreds of gigabytes.

Four techniques bring that down, and real assemblers combine them. First, partition k-mers by their minimizer, the smallest m-mer they contain, so consecutive k-mers usually land in the same bucket and each bucket can be counted on disk and in memory separately. KMC works this way. Second, filter singletons before inserting them, often with a Bloom filter that admits a k-mer to the real table only on its second sighting. Third, represent the graph probabilistically: Minia stores solid k-mers in a Bloom filter and keeps a small exact table of the false positives that would create wrong branches. Fourth, use a succinct representation: MEGAHIT's BOSS structure encodes the graph in a few bits per edge with rank and select queries. The underlying sketches are explained in Bloom filter variants and the Count-Min sketch.

Failure modes

Adapter and quality residue. Untrimmed adapters create high-count k-mers that link unrelated contigs. Trim adapters and low-quality tails before counting, and check that the histogram's error spike shrinks.

Threshold set on the label, not the histogram. A fixed count threshold that suits 60x data deletes real sequence at 15x. Read the trough from your own histogram every time.

Heterozygosity mistaken for errors. Highly heterozygous diploid genomes produce a half-coverage peak and dense bubbles. Popping all of them yields a mosaic of the two haplotypes, and leaving them yields a duplicated assembly that is almost twice the true size. Use a haplotype-aware tool when the half-coverage peak is large.

Strand bugs in home-grown code. Forgetting to add both orientations of a canonical k-mer as edges produces a graph where contigs stop at every strand change. Test on a genome whose reads come from both strands and check that the contig matches the genome or its reverse complement.

Trade-offs

When choosing a string index for searches against the finished contigs, rather than for building them, suffix arrays are the usual tool; see suffix arrays in depth.

What to do next

1. Run the toy assembler above on the example genome, then change the error position and confirm that a middle error produces a bubble rather than a tip.
2. Count k-mers on a real sample with a dedicated counter and plot the histogram before choosing any threshold.
3. Compute k-mer coverage and the error-free fraction for your read length and three candidate values of k.
4. Assemble with a multi-k assembler, then report contig N50, total length and the fraction of reads that map back.
5. Inspect the assembly graph in a viewer and find the largest unresolved repeat.
6. If contiguity matters, add paired-end or long-read evidence rather than tuning k further.