Genome assembly

De novo assembly reconstructs a genome from reads without a reference: it is the jigsaw puzzle at the heart of genomics. Because every read is a short window onto a long molecule, assembly works by finding overlaps between reads and following them into longer contiguous sequences, or contigs. This lesson explains the two main families of assembly algorithms, the overlap-layout-consensus approach and de Bruijn graphs, why repeats are the fundamental obstacle, how assemblies are scaffolded and evaluated, and how long reads have made complete, haplotype-resolved genomes routine.

The assembly problem

In its simplest formulation, assembly asks for the shortest string that contains every read as a substring, the shortest common superstring. That problem is NP-hard, which means that no efficient algorithm is known for the general case, and it is also the wrong model: repeats in real genomes mean that the shortest superstring collapses repeated copies into one. Practical assemblers instead build a graph of how reads or k-mers fit together and then simplify it, using coverage and read connections to decide how to traverse it.

Repeats are the central difficulty. If a sequence of length \(R\) occurs twice in a genome, and every read is shorter than \(R\), then no read connects a unique flanking region on one side of a repeat copy to the unique region on the other side, and the assembler cannot tell which left flank goes with which right flank. The graph contains the information that the repeat exists and how many copies there may be, but not how the copies are arranged. The assembly therefore breaks at the repeat, or worse, joins the wrong flanks in a misassembly. The longest common repeats set the read length needed: bacterial rRNA operons of 5 to 7 kb, insertion sequences of 1 to 2 kb, human LINE-1 elements of about 6 kb, segmental duplications of tens to hundreds of kilobases and centromeric arrays of megabases.

Overlap–layout–consensus

The overlap-layout-consensus (OLC) approach, used in the Sanger era by the Celera assembler, computes overlaps between all pairs of reads, builds an overlap graph with reads as nodes and overlaps as edges, finds a layout, an ordering of reads that is a path through the graph, and computes a consensus sequence from the aligned reads. Finding a path that visits every node exactly once, a Hamiltonian path, is computationally hard in general, but assemblers simplify the graph first: they remove transitively inferable edges (if A overlaps B and B overlaps C, the edge from A to C is redundant) to produce a string graph, and they remove reads contained within others. The overlap step, all-against-all comparison, was once prohibitive for millions of reads, but fast methods based on minimizers make it practical for long reads.

OLC suits long reads, where reads are few enough to compare and long enough that overlaps are informative. Modern long-read assemblers descend from it: Canu, for noisy long reads; hifiasm, which builds string graphs from accurate HiFi reads and keeps the two haplotypes of a diploid apart; and Flye, which first builds a repeat graph from approximate overlaps and then resolves repeats with the reads that span them. With HiFi reads of 15 to 25 kb and ultra-long nanopore reads, assemblers now routinely produce chromosome-scale contigs for many species.

De Bruijn graph assembly

When short-read sequencers produced hundreds of millions of reads, all-against-all overlap became impractical, and a different idea, proposed for sequencing by Pavel Pevzner and colleagues in 2001, took over. Break every read into k-mers. Build a graph in which each distinct \((k-1)\)-mer is a node and each k-mer is an edge from its prefix \((k-1)\)-mer to its suffix \((k-1)\)-mer. Consecutive k-mers in the genome share \(k - 2\) bases and so form a walk through consecutive edges, and the genome itself corresponds to a path that uses every edge once, an Eulerian path. Unlike the Hamiltonian path, an Eulerian path can be found in time proportional to the number of edges, an observation going back to Leonhard Euler's 1736 solution of the Königsberg bridges problem.

Real de Bruijn graphs are messy. A sequencing error creates up to \(k\) erroneous k-mers that form a short dead-end branch, a tip, or a short parallel path, a bubble, rejoining the main path; heterozygous sites in a diploid also create bubbles. Assemblers such as Velvet, SPAdes and MEGAHIT remove tips and pop bubbles using coverage: error branches have low coverage. Repeats longer than \(k\) collapse into a single path with entry and exit branches. Unambiguous stretches of the simplified graph, where nodes have one way in and one way out, become contigs, called unitigs. The choice of \(k\) balances repeat resolution (larger \(k\) resolves more repeats) against connectivity (larger \(k\) means fewer error-free k-mers and more gaps at low coverage), so SPAdes builds graphs with several values of \(k\) and combines them.

# Build a de Bruijn graph from error-free reads of a small random genome and find an Eulerian path.
import random
from collections import defaultdict, Counter
random.seed(41)
genome = ''.join(random.choice('ACGT') for _ in range(400))
reads = [genome[i:i + 50] for i in range(0, len(genome) - 49, 7)] + [genome[-50:]]
def assemble(reads, k):
    kmers = set()
    for r in reads:
        for i in range(len(r) - k + 1):
            kmers.add(r[i:i + k])
    graph = defaultdict(list)
    indeg, outdeg = Counter(), Counter()
    for km in kmers:
        a, b = km[:-1], km[1:]
        graph[a].append(b)
        outdeg[a] += 1
        indeg[b] += 1
    starts = [n for n in graph if outdeg[n] - indeg[n] == 1]
    start = starts[0] if starts else next(iter(graph))
    stack, path = [start], []          # Hierholzer's algorithm
    while stack:
        v = stack[-1]
        if graph[v]:
            stack.append(graph[v].pop())
        else:
            path.append(stack.pop())
    path.reverse()
    return path[0] + ''.join(n[-1] for n in path[1:]), len(kmers), len(starts)

for k in (6, 9, 15):
    seq, n_edges, n_starts = assemble(reads, k)
    print(f"k={k:2d}: {n_edges:3d} distinct k-mers, path length {len(seq)}, matches genome: {seq == genome}")

k= 6: 385 distinct k-mers, path length 372, matches genome: False
k= 9: 392 distinct k-mers, path length 400, matches genome: True
k=15: 386 distinct k-mers, path length 400, matches genome: True

With \(k = 6\), only 4,096 possible 6-mers exist and a 400-base random sequence repeats some of its 5-mers by chance, so repeated 6-mers collapse into single edges and the graph no longer records how many times each was used; the path found is 28 bases short and rearranged relative to the genome. Once \(k\) is large enough that every \((k-1)\)-mer is unique, the Eulerian path is unique and reconstructs the genome exactly. The same logic applies to real genomes, except that their repeats are not chance coincidences of short words but genuine copies thousands of bases long, which no practical \(k\) can resolve with short reads.

Scaffolding and gap closing

Contigs end at repeats, coverage gaps and unresolved heterozygosity. Scaffolding orders and orients contigs and estimates the gaps between them, using information that spans longer distances than the reads used for contig building. Paired-end reads with inserts of a few hundred bases, and formerly mate-pair libraries with inserts of several kilobases, link contigs whose ends they span. Long reads link contigs across repeats. Chromosome conformation capture (Hi-C) exploits the fact that DNA segments on the same chromosome contact each other far more often than segments on different chromosomes, with contact frequency falling with distance; mapping Hi-C read pairs to contigs reveals which contigs belong to the same chromosome and in what order. Optical maps, images of long DNA molecules labeled at specific sequence motifs, provide independent long-range structure. Gaps of estimated length are filled with runs of N in the scaffold sequence.

Evaluating an assembly

How good is an assembly? Contiguity is measured by the N50: sort contigs from longest to shortest and add their lengths until the running total reaches half the total assembly size; the length of the contig at which this happens is the N50, and the number of contigs needed is the L50. Because N50 depends on the assembly's own total size, which may be inflated by duplicated haplotypes or contamination, the NG50 uses half the estimated genome size instead. Contiguity is not correctness: an assembler that joins contigs aggressively achieves a high N50 with misassemblies.

# Compute N50, L50 and NG50 for an invented set of contig lengths (bp).
contigs = [2_450_000, 1_310_000, 980_000, 640_000, 410_000, 220_000, 150_000, 95_000,
           60_000, 41_000, 30_000, 22_000, 15_000, 9_000, 6_500, 4_000, 2_800, 2_100]
def nx(lengths, target):
    total = 0
    for i, n in enumerate(sorted(lengths, reverse=True), start=1):
        total += n
        if total >= target:
            return n, i
    return None, None
size = sum(contigs)
n50, l50 = nx(contigs, size / 2)
ng50, lg50 = nx(contigs, 7_200_000 / 2)       # estimated genome size from k-mers: 7.2 Mb
print(f"total assembly size {size:,} bp in {len(contigs)} contigs")
print(f"N50 {n50:,} bp (L50 = {l50}),  NG50 {ng50:,} bp (LG50 = {lg50})")
print(f"assembly covers {size / 7_200_000:.1%} of the estimated genome size")

total assembly size 6,447,400 bp in 18 contigs
N50 1,310,000 bp (L50 = 2),  NG50 1,310,000 bp (LG50 = 2)
assembly covers 89.5% of the estimated genome size

Completeness and correctness need other measures. BUSCO searches the assembly for a set of genes expected to be present in single copy in nearly all members of a lineage, such as several thousand conserved genes for mammals; a good assembly finds above 95 percent of them complete, few fragmented, and few duplicated (many duplicates suggest uncollapsed haplotypes). k-mer-based tools such as Merqury compare k-mers in the assembly with those in accurate reads: k-mers found in the assembly but not in reads indicate base errors, giving a consensus quality value (QV) on the Phred scale, where QV 40 means one error per 10,000 bases and QV 60 one per million. Mapping reads back reveals collapsed repeats (twice the expected coverage) and misjoins, and comparison with a related genome or with Hi-C contact maps reveals large misassemblies.

Haplotypes, polishing and contamination

A diploid genome has two haplotypes, and older assemblers produced a single mosaic sequence that switched between them. Haplotype-resolved assembly keeps them apart. Trio binning sequences both parents and sorts the child's long reads by k-mers unique to the mother or father before assembling each set separately. Without parents, hifiasm uses the accuracy of HiFi reads to separate haplotypes within the graph, and Hi-C or long-range phasing assigns blocks to chromosomes. Polishing corrects residual base errors by aligning reads to the assembly and taking a consensus, important for noisier reads. Finally, contamination must be removed: contigs from bacteria in a tissue sample, from the sequencing control PhiX, or from adapter sequence. Submitting a genome to a public database now includes automated contamination screening.