bio-applied-genome-assembly

name	bio-applied-genome-assembly
description	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.
tool_type	python
primary_tool	Python

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.

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 = (overlaps, key= x: x[])
        merged = contigs[i] + contigs[j][ov:]
        contigs = [c  idx, c  (contigs)  idx   (i, j)] + [merged]
     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  edges:
            in_deg[dest] += 

    all_nodes = (out_deg) | (in_deg)
    start = ((n  n  all_nodes  out_deg.get(n, ) - in_deg.get(n, ) == ), )
     start  :
        start = ((out_deg))

    stack, path = [start], []
     stack:
        v = stack[-]
         adj.get(v):
            stack.append(adj[v].pop())
        :
            path.append(stack.pop())
    path.reverse()
     path


 () -> :
    
      path:
         
     path[] + .join(node[-]  node  path[:])

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

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

# 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