# Bio Applied Genome Assembly

> Implement OLC and de Bruijn assembly algorithms, compute N50/L50/NG50 stats, and run SPAdes/Flye/hifiasm on Illumina/HiFi/ONT reads. Use for k-mer graphs, comparing assemblers, or a FASTQ-to-contigs pipeline.

- Skill: `pavel-kravchenko/bio-applied-genome-assembly` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-genome-assembly`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-genome-assembly/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: DevOps & Infra
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

| Algorithm | Reads | Key idea | Complexity |
|---|---|---|---|
| OLC (Overlap-Layout-Consensus) | Long reads (PacBio, ONT, Sanger) | Hamiltonian path in overlap graph | O(n²) overlaps |
| De Bruijn graph | Short reads (Illumina) | Eulerian path in k-mer graph | O(n·L) |

**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.

```python
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
```

```python
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).

```python
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

| Assembler | Read type | Algorithm | Best for |
|---|---|---|---|
| **SPAdes** | Illumina | Multi-k de Bruijn | Bacteria, small eukaryotes, metagenomes |
| **Canu** | PacBio CLR / ONT | OLC + correction | Large genomes, high repeat content |
| **Flye** | PacBio CLR / ONT / HiFi | Repeat graph | Fast, tolerates high error rates |
| **hifiasm** | PacBio HiFi | Graph-based | Haplotype-resolved assembly |
| **verkko** | HiFi + ONT ultra-long | Graph-based | Telomere-to-telomere assemblies |

```bash
# 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