Bio Applied Advanced Ngs

Genome Assembly and Advanced NGS

When to Use

• Choosing between reference mapping and de novo assembly for a new organism or structural-variant discovery.
• Picking an assembler (SPAdes, Flye, hifiasm, Canu, verkko) based on read type (Illumina, ONT, PacBio CLR/HiFi).
• Explaining or debugging why a de Bruijn / OLC assembly graph broke (tips, bubbles, repeat tangles).
• Choosing k-mer size for a de Bruijn assembler or diagnosing k-mer-related fragmentation.
• Computing/comparing assembly contiguity (N50, L50, NG50, N90) across assemblies or genome-size assumptions.

Version Compatibility

SPAdes ≥4.0, Flye ≥2.9, hifiasm ≥0.19, QUAST ≥5.2, BUSCO ≥5.5, Python ≥3.10 (uses list[str]/dict[str, list[str]] type hints).

Prerequisites

• Bioconda packages: spades, flye, hifiasm, quast, busco.
• Python stdlib only for the algorithmic parts (collections.defaultdict).
• Concepts: FASTQ/FASTA basics, sequencing coverage, k-mers. See bio-sequence-io-read-sequences for I/O.

OLC Assembly (long reads)

Goal: reconstruct a sequence from overlapping reads without a reference. Approach: find all suffix–prefix overlaps, then greedily merge the pair with the longest overlap until none remain. O(n²) per round — only tractable for hundreds to thousands of long reads, not billions of short reads.

def find_overlaps(reads: list[str], min_overlap: int = 4) -> list[tuple]:
    """Find all suffix-prefix overlaps between reads (i's suffix == j's prefix)."""
    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


reads = ["AGCTACAGTATGCT", "TACAGTATGCTTAT", "GTATGCTTATCTGA",
         "TGCTTATCTGATAC", "TGATACCTTAGCCA"]
result = greedy_assemble_olc(reads, min_overlap=5)
assert result[0] == "AGCTACAGTATGCTTATCTGATACCTTAGCCA"

De Bruijn Graph Assembly (short reads)

Goal: reconstruct a sequence from millions of short reads in near-linear time. Approach: decompose reads into overlapping k-mers, build a graph where nodes are (k-1)-mers and edges are k-mers, then walk an Eulerian path (Hierholzer's algorithm) through the graph. Sequencing errors create dead-end "tips"; heterozygous SNPs create "bubbles"; repeats longer than k create branch points that need paired-end or long-read scaffolding to resolve.

from collections import defaultdict

def build_debruijn_graph(reads: list[str], k: int) -> dict[str, list[str]]:
    """Build a de Bruijn graph. Nodes: (k-1)-mers. Edges: k-mers."""
    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]:
    """Find an Eulerian path via Hierholzer's algorithm.

    Starts at the node with out_degree - in_degree == 1 (path start);
    falls back to any node with edges for an Eulerian circuit.
    """
    adj = {node: list(edges) for node, edges in graph.items()}
    out_deg = {node: len(edges) for node, edges in adj.items()}
    in_deg = 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), 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 from a de Bruijn node path."""
    if not path:
        return ""
    return path[0] + "".join(node[-1] for node in path[1:])


target = "AGCTACAGTATGCTTATCTGATAC"
reads = [target[i:i + 10] for i in range(0, len(target) - 10 + 1, 2)]
graph = build_debruijn_graph(reads, k=9)
assembled = path_to_sequence(eulerian_path(graph))
assert assembled == target

Assembly QC: N50 / L50 / NG50

Goal: quantify assembly contiguity to compare tools/parameters objectively. Approach: sort contig lengths descending; N50 is the length at which the cumulative sum crosses half the total assembly length (or half the known genome size for NG50); L50 is the number of contigs needed to reach that point.

def calculate_n50_l50(contig_lengths: list[int]) -> tuple[int, int]:
    """Calculate N50 and L50 from a list of contig lengths."""
    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:
    """NG50: like N50, but the threshold is half the estimated genome size."""
    cumsum = 0
    for length in sorted(contig_lengths, reverse=True):
        cumsum += length
        if cumsum >= genome_size / 2:
            return length
    return 0  # assembly total shorter than the genome-size estimate

CLI Usage

# 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 (ONT raw, PacBio CLR, or HiFi via --pacbio-hifi)
flye --nano-raw reads.fastq.gz --genome-size 5m --out-dir flye_out/ --threads 8

# hifiasm -- PacBio HiFi, haplotype-resolved
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

Assembler Reference

SPAdes improvements over a naive de Bruijn graph: runs multiple k values (e.g. 21, 33, 55, 77) and merges graphs; discards low-count k-mers as errors; clips tips (dead-end branches); pops bubbles (merges heterozygous SNP paths); scaffolds contigs using paired-end insert-size distributions.

Pitfalls

• k-mer size is critical: too small → ambiguous graph, unresolvable repeats; too large → breaks at sequencing errors. Practical range k=31–127 for Illumina; SPAdes uses multiple k simultaneously.
• OLC is O(n²) — not for short reads: pairwise overlap computation is prohibitive for billions of Illumina reads; use de Bruijn graphs there instead.
• Repeat content breaks assemblies: repeats longer than the read length create irresolvable branches. Long reads (ONT/HiFi) are essential for repeat-rich genomes.
• N50 is not correctness: a high N50 can result from chimeric contigs. Always run QUAST/BUSCO, not just N50, before trusting an assembly.
• Coverage floors: de novo assembly typically needs ≥50x coverage; below this, contigs break at low-coverage regions.
• hifiasm GFA output: hifiasm emits GFA, not FASTA directly — convert with awk (see CLI example) before downstream tools that expect FASTA.

See Also

• bio-genome-assembly-short-read-assembly — SPAdes/Illumina workflow details.
• bio-genome-assembly-long-read-assembly — Flye/Canu ONT workflow details.
• bio-genome-assembly-hifi-assembly — hifiasm/verkko HiFi workflow details.
• bio-genome-assembly-assembly-qc — QUAST/BUSCO evaluation in depth.