bio-genome-assembly-short-read-assembly

What this skill does

## Version Compatibility

Reference examples tested with: SPAdes 4.0+, MEGAHIT 1.2+, Unicycler 0.5+, ABySS 2.3+, GenomeScope2 2.0+, KMC 3.2+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags

SPAdes 4.0 (June 2024) requires Python >=3.8, defaults to GFA v1.2 output with circular-path tags (`TP:Z:circular`), and is stated by the maintainers to be the last major feature release - pin the version in methods, because k auto-selection and the `--isolate`/`--careful` recommendation flipped across 3.x. MEGAHIT preset k-lists (`meta-sensitive`, `meta-large`) have changed across versions - confirm against `megahit -h` before hard-coding. If code throws an error, introspect the installed tool and adapt rather than retrying.

# Short-Read Assembly

**"Assemble a genome from Illumina reads"** -> Build a multi-k de Bruijn graph from short reads and walk it into contigs, knowing the contiguity ceiling is set by the genome's repeat structure, not the depth.
- CLI: `spades.py --isolate -1 R1.fq.gz -2 R2.fq.gz -o out -t 16` (bacterial isolate, the genuine sweet spot)
- CLI: `megahit -1 R1.fq.gz -2 R2.fq.gz -o out -t 16` (huge/low-memory, metagenome)
- CLI: profile first (see below): `genomescope2 -i kmer.hist -o gs -k 21 -p 2` (size/het/repeat before parameters)

## The Single Most Important Modern Insight -- Short Reads Cannot Span Repeats, So Contiguity Is Capped by the Genome, Not the Depth

A short read (~150 bp) or paired insert (~300-600 bp) cannot resolve any repeat longer than itself: the de Bruijn graph enters a repeat from multiple unique flanks, traverses an identical internal path, and the assembler's only safe move is to break the contig at the repeat boundary or collapse the copies. Genomic repeats - rRNA operons (~5 kb), transposons, segmental duplications, satellites - are routinely kilobases to megabases. So a short-read assembly is **structurally fragmented at every long repeat**, and that is biology, not a software defect. Three load-bearing corollaries:

1. **N50 plateaus at the repeat structure, not the sequencing depth.** Past ~50-100x for an isolate (~50-60x for a eukaryote), more Illumina does NOT lengthen contigs - it only adds error-correction signal, and above ~150x it actively hurts by amplifying GC bias and duplicate-driven coverage distortion. The reflex "add more reads to get a better assembly" is wrong: the limit is read length.

2. **The assembler is almost never the bottleneck - the input DNA and the genome's repeat/heterozygosity structure are.** A clean, high-molecular-weight, PCR-free library at the right depth assembles beautifully with defaults; a degraded/contaminated library assembles badly with every setting and every k. The expert's first move on a bad assembly is to look at the k-mer spectrum, GC-vs-coverage, and duplication rate (GenomeScope2 profiling, see the profiling section below), not to swap assemblers or sweep k.

3. **Short reads are no longer the frontier instrument for a finished genome.** For any genome that must be complete or finished, the answer is long reads that span the repeats (`-> long-read-assembly`, `-> hifi-assembly`), with short reads relegated to polishing, k-mer-spectrum QC, and the cheap bacterial-isolate/surveillance workhorse. A skill that lets an agent attempt a de novo short-read assembly of a large, repeat-rich, heterozygous eukaryote and then blame parameters is doing harm.

## Tool Taxonomy

| Tool | Citation | Role | When |
|------|----------|------|------|
| SPAdes | Bankevich 2012 *J Comput Biol* | multi-k de Bruijn assembler with mode dispatch | bacterial/fungal/small-euk isolate; the field default |
| MEGAHIT | Li 2015 *Bioinformatics* | succinct (compressed) de Bruijn, ultra-low memory | huge short-read datasets, metagenomes, memory-constrained |
| Unicycler | Wick 2017 *PLoS Comput Biol* | SPAdes wrapper + graph bridging + circularization | bacterial finishing; hybrid Illumina+long-read |
| MaSuRCA | Zimin 2013 *Bioinformatics* | super-reads + mega-reads (hybrid OLC/DBG) | large eukaryotic genomes with mixed data |
| ABySS 2.0 | Jackman 2017 *Genome Res* | Bloom-filter de Bruijn, MPI-parallel | large genomes on memory-constrained HPC |
| Platanus / Platanus-allee | Kajitani 2014 *Genome Res* | bubble-aware DBG for high heterozygosity | highly heterozygous diploids (a real niche) |
| GenomeScope2 / Smudgeplot | Ranallo-Benavidez 2020 *Nat Commun* | k-mer-histogram model: size/het/repeat/ploidy | run FIRST; sets the ceiling reference-free |
| Velvet / SOAPdenovo2 | Zerbino 2008 *Genome Res* / Luo 2012 *GigaScience* | single-k DBG, superseded | reproduce old papers only; do not start new projects |

Every mainstream short-read assembler is a de Bruijn assembler because all-vs-all overlap (OLC) is infeasible for hundreds of millions of short reads. The single most consequential parameter is k: small k over-connects (collapses repeats, chimeras), large k is more unique (resolves repeats up to length k) but demands higher coverage and is error-fragile. No single k is optimal everywhere in a genome, which is exactly why SPAdes and MEGAHIT iterate over a k-series - do NOT hand-pick a single k.

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Reads not yet profiled | GenomeScope2 on a k-mer histogram (profiling section below) | the spectrum sets size/het/repeat ceiling before any parameter |
| Reads not yet QC'd | -> read-qc/quality-reports, read-qc/adapter-trimming | over-trimming creates coverage holes that fragment the graph |
| Bacterial/archaeal/viral isolate, Illumina only | SPAdes `--isolate` | multi-k DBG + paired-end repeat resolution; the genuine sweet spot |
| Small bacterial genome, want mismatch/indel cleanup | SPAdes `--careful` (small genomes ONLY) | runs MismatchCorrector; NOT for large eukaryotes; incompatible with --meta/--rna |
| Bacterial isolate, want a finished/circular assembly | Unicycler (short-only or hybrid) | graph bridging + circularization + rotation to dnaA |
| Single-cell / MDA-amplified | SPAdes `--sc` | designed for the wildly uneven MDA coverage |
| Highly heterozygous diploid, short reads | Platanus-allee, then purge | bubble-aware; but consider -> hifi-assembly instead |
| Metagenome (recover MAGs) | -> metagenome-assembly (metaSPAdes / MEGAHIT) | community co-assembly; MAG recovery, not single-genome N50 |
| Huge dataset, RAM-constrained | MEGAHIT | succinct DBG, tiny memory footprint |
| Large/repeat-rich/heterozygous eukaryote, finished genome wanted | -> hifi-assembly or -> long-read-assembly | short reads mathematically cannot span the repeats |
| Assembly built, judging quality | -> assembly-qc (NG50 + auN + BUSCO + Merqury QV) | never report N50 alone; contiguity is not correctness |

## SPAdes -- the dominant short-read assembler (the mode is a graph choice, not a cosmetic flag)

```bash
spades.py --isolate -1 R1.fq.gz -2 R2.fq.gz -o out -t 16 -m 64   # recommended default for isolates; does NOT run MismatchCorrector
spades.py --careful -1 R1.fq.gz -2 R2.fq.gz -o out               # SMALL genomes only; runs MismatchCorrector; NOT eukaryotes; not with --meta/--rna
spades.py --sc      -1 R1.fq.gz -2 R2.fq.gz -o out               # single-cell/MDA (default k 21,33,55)
spades.py --plasmid -1 R1.fq.gz -2 R2.fq.gz -o out              # plasmidSPAdes (coverage-based extraction)
spades.py -1 R1.fq.gz -2 R2.fq.gz -o out -k 21,33,55,77         # explicit multi-k LIST (odd, <= read length) only if overriding
spades.py --only-assembler -1 R1.fq.gz -2 R2.fq.gz -o out       # skip BayesHammer (reads pre-corrected, or het/pooled data)
```

Picking the wrong mode is the most common SPAdes error. `--isolate` and `--careful` solve different problems: modern SPAdes pushes `--isolate` as the isolate default (tuned for high, even coverage), while `--careful` runs BWA-based MismatchCorrector to cut mismatches/short indels