| name | bio-applied-advanced-ngs |
| description | 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. |
| tool_type | bash |
| primary_tool | SPAdes |
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[])
merged = contigs[i] + contigs[j][ov:]
contigs = [c idx, c (contigs) idx (i, j)] + [merged]
contigs
reads = [, , ,
, ]
result = greedy_assemble_olc(reads, min_overlap=)
result[] ==
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 = (out_deg) | (in_deg)
start = ((n n all_nodes
out_deg.get(n, ) - in_deg.get(n, ) == ), ((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[:])
target =
reads = [target[i:i + ] i (, (target) - + , )]
graph = build_debruijn_graph(reads, k=)
assembled = path_to_sequence(eulerian_path(graph))
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
CLI Usage
spades.py -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -o spades_out/ -t 8
spades.py --meta -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -o metaspades_out/
flye --nano-raw reads.fastq.gz --genome-size 5m --out-dir flye_out/ --threads 8
hifiasm -o assembly -t 16 hifi_reads.fastq.gz
awk '/^S/{print ">"$2"\n"$3}' assembly.bp.p_ctg.gfa > assembly.p_ctg.fasta
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.
