# Bio Applied Advanced Ngs

> Assemble genomes de novo: greedy OLC, de Bruijn graph/Eulerian path, N50/L50/NG50 stats, SPAdes/Flye/hifiasm CLI usage. Use when choosing k-mer size, picking an assembler for Illumina/ONT/HiFi reads, or scoring contiguity.

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

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

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

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

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

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