Bio Applied Genome Assembly

Genome Assembly

When to Use

• Explaining or implementing how genome assemblers work (OLC vs. de Bruijn graph, Eulerian path)
• Choosing k-mer size or an assembler (SPAdes, Canu, Flye, hifiasm, verkko) for a given read type
• Computing/interpreting assembly-quality metrics (N50, L50, NG50) to compare two assemblies
• Building a FASTQ → contigs de novo assembly pipeline for bacteria, eukaryotes, or metagenomes
• Diagnosing why an assembly is fragmented (low coverage, repeats, wrong k)

Version Compatibility

SPAdes ≥3.15, Flye ≥2.9, hifiasm ≥0.19, Canu ≥2.2, QUAST ≥5.2, BUSCO ≥5.4. Python ≥3.10 (uses list[str] / X | None type hints). No third-party Python packages needed for the algorithm code below — pure stdlib (collections.defaultdict).

Prerequisites

• Command-line tools on PATH: spades.py, flye, hifiasm, quast.py, busco (install via conda install -c bioconda spades flye hifiasm quast busco)
• Concepts: k-mers, graph traversal (Eulerian/Hamiltonian paths), FASTQ/FASTA formats
• Related upstream skill: bio-read-qc-fastp-workflow for trimming/QC before assembly

Algorithm Comparison

Goal: reconstruct a genome sequence from a set of overlapping reads. Approach: for long reads, find pairwise suffix-prefix overlaps and greedily merge the pair with the longest overlap (OLC). For short reads, this is O(n²) and infeasible — instead break reads into k-mers, build a de Bruijn graph where nodes are (k-1)-mers and edges are k-mers, then walk an Eulerian path through it.

from collections import defaultdict


def find_overlaps(reads: list[str], min_overlap: int = 4) -> list[tuple]:
    """Find all suffix-prefix overlaps between reads (all-vs-all, O(n^2))."""
    overlaps = []
    for i, r1 in enumerate(reads):
        for j, r2 in enumerate(reads):
            if i == j:
                continue
            max_overlap = min(len(r1), len(r2))
            for k in range(max_overlap, min_overlap - 1, -1):
                if r1[-k:] == r2[:k]:
                    overlaps.append((i, j, k))
                    break
    return overlaps


def greedy_assemble_olc(reads: list[str], min_overlap: int = 4) -> list[str]:
    """Greedy OLC assembler: repeatedly merge the pair with the longest overlap."""
    contigs = list(reads)
    while True:
        overlaps = find_overlaps(contigs, min_overlap)
        if not overlaps:
            break
        i, j, ov = max(overlaps, key=lambda x: x[2])
        merged = contigs[i] + contigs[j][ov:]
        contigs = [c for idx, c in enumerate(contigs) if idx not in (i, j)] + [merged]
    return contigs

from collections import defaultdict


def build_debruijn_graph(reads: list[str], k: int) -> dict[str, list[str]]:
    """Build a de Bruijn graph. Nodes are (k-1)-mers, edges are k-mers.

    Returns adjacency list: {left_kmer: [right_kmer, ...]}
    """
    graph = defaultdict(list)
    for read in reads:
        for i in range(len(read) - k + 1):
            kmer = read[i:i + k]
            graph[kmer[:-1]].append(kmer[1:])
    return dict(graph)


def eulerian_path(graph: dict[str, list[str]]) -> list[str] | None:
    """Find an Eulerian path in a directed graph using Hierholzer's algorithm."""
    adj = {node: list(edges) for node, edges in graph.items()}
    out_deg = {node: len(edges) for node, edges in adj.items()}
    in_deg: dict[str, int] = defaultdict(int)
    for edges in adj.values():
        for dest in edges:
            in_deg[dest] += 1

    all_nodes = set(out_deg) | set(in_deg)
    start = next((n for n in all_nodes if out_deg.get(n, 0) - in_deg.get(n, 0) == 1), None)
    if start is None:
        start = next(iter(out_deg))

    stack, path = [start], []
    while stack:
        v = stack[-1]
        if adj.get(v):
            stack.append(adj[v].pop())
        else:
            path.append(stack.pop())
    path.reverse()
    return path


def path_to_sequence(path: list[str]) -> str:
    """Reconstruct the assembled sequence by concatenating overlapping (k-1)-mers."""
    if not path:
        return ""
    return path[0] + "".join(node[-1] for node in path[1:])

de Bruijn graph complications in practice:

• Sequencing errors → spurious k-mers → "tip" branches (dead ends)
• Diploid SNPs → "bubbles" (two paths with same entry/exit)
• Repeats → complex nodes with multiple in/out edges

Assembly Quality Metrics

Goal: quantify and compare assembly contiguity so two assemblies (or assembler settings) can be judged objectively. Approach: sort contig lengths descending and find the length at which cumulative sum crosses 50% (N50/L50) of the assembly total, or 50% of the estimated genome size (NG50, more comparable across differently-sized assemblies of the same genome).

def calculate_n50_l50(contig_lengths: list[int]) -> tuple[int, int]:
    """Return (N50, L50): N50 is the contig length at which cumulative
    length crosses half the assembly total; L50 is the number of contigs
    needed to reach it (fewer = more contiguous).
    """
    lengths = sorted(contig_lengths, reverse=True)
    threshold = sum(lengths) / 2
    cumsum = 0
    for i, length in enumerate(lengths):
        cumsum += length
        if cumsum >= threshold:
            return length, i + 1
    return lengths[-1], len(lengths)


def calculate_ng50(contig_lengths: list[int], genome_size: int) -> int:
    """N50 using the estimated genome size (not assembly total) as denominator."""
    lengths = sorted(contig_lengths, reverse=True)
    threshold = genome_size / 2
    cumsum = 0
    for length in lengths:
        cumsum += length
        if cumsum >= threshold:
            return length
    return 0  # assembly total < genome_size estimate

Assembler Reference

# SPAdes -- Illumina paired-end
spades.py -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -o spades_out/ -t 8

# SPAdes -- metagenome
spades.py --meta -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -o metaspades_out/

# Flye -- long reads
flye --nano-raw reads.fastq.gz --genome-size 5m --out-dir flye_out/ --threads 8

# hifiasm -- PacBio HiFi (outputs assembly.bp.p_ctg.gfa; convert to FASTA with awk)
hifiasm -o assembly -t 16 hifi_reads.fastq.gz
awk '/^S/{print ">"$2"\n"$3}' assembly.bp.p_ctg.gfa > assembly.p_ctg.fasta

# Assembly QC
quast.py contigs.fasta -r reference.fasta -o quast_out/
busco -i contigs.fasta -l bacteria_odb10 -o busco_out/ -m genome

SPAdes Graph Simplification Steps

1. Runs multiple k values (e.g., 21, 33, 55, 77) and merges graphs
2. Error correction: k-mers below coverage threshold are discarded
3. Tip clipping: dead-end branches shorter than threshold
4. Bubble popping: merges alternative paths caused by SNPs
5. Paired-end scaffolding: links contigs using insert size information

Pitfalls

• k-mer size is critical: Too small (k=2) → ambiguous graph, can't resolve repeats. Too large → breaks at sequencing errors, poor connectivity. Practical: k=31–127 for Illumina; SPAdes uses multiple k simultaneously.
• Repeat content breaks assemblies: Repeats longer than the read length create irresolvable branches in the de Bruijn graph. Long reads (PacBio/ONT) are essential for repeat-rich genomes.
• Short reads cannot use OLC: O(n²) pairwise overlaps are computationally prohibitive for billions of Illumina reads — use de Bruijn graph instead.
• N50 is not a measure of correctness: High N50 can result from chimeric contigs. Always run BUSCO or QUAST/CheckM after assembly.
• N50 vs NG50: N50 uses assembly total length as the denominator; NG50 uses the estimated genome size. If the assembly is incomplete, N50 will look better than NG50 — always report NG50 when comparing assemblies against a known genome size.
• Coverage matters: De novo assembly requires ~50–100x coverage. Lower coverage creates gaps and breaks contigs at low-coverage regions.

See Also

• bio-genome-assembly-short-read-assembly — production SPAdes workflows for Illumina data
• bio-genome-assembly-hifi-assembly — hifiasm/verkko workflows for PacBio HiFi
• bio-genome-assembly-assembly-qc — BUSCO/QUAST QC pipelines in depth
• bio-workflows-genome-assembly-pipeline — end-to-end FASTQ-to-assembly pipeline