Bio Genome Assembly Short Read Assembly

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

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

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

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 on small genomes only - it is explicitly slow, memory-heavy, and low-benefit on medium/large eukaryotes, and is incompatible with --meta and --rna. --isolate and --careful are mutually exclusive - SPAdes aborts (cannot specify --careful in isolate mode) if both are given; pick one. --meta, --rna, --bio dispatch to genuinely different pipelines (metaSPAdes, rnaSPAdes, biosyntheticSPAdes) with different graph models - they are NOT genome assembly. BayesHammer (the built-in Illumina error-corrector) runs by default; --only-assembler skips it. Outputs: contigs.fasta, scaffolds.fasta, assembly_graph.gfa. Always report contig metrics, not just scaffold metrics (scaffold gaps are estimated N-runs, not sequence).

MEGAHIT, Unicycler, and the others

megahit -1 R1.fq.gz -2 R2.fq.gz -o out -t 16                    # default k 21,29,39,59,79,99,119,141; min-count 2 (drops singleton error k-mers)
megahit -1 R1.fq.gz -2 R2.fq.gz --presets meta-sensitive -o out # min-count 1; verify k-list against megahit -h
unicycler -1 R1.fq.gz -2 R2.fq.gz -o out -t 16                  # short-read-only bacterial finishing
unicycler -1 R1.fq.gz -2 R2.fq.gz -l long.fq.gz -o out          # hybrid (short for accuracy + long to bridge repeats)
abyss-pe name=asm k=96 B=2G in='R1.fq.gz R2.fq.gz'              # B = Bloom-filter size (the 2.0 low-memory mode); single k

MEGAHIT's succinct DBG gives a tiny memory footprint at some cost in contiguity/accuracy versus metaSPAdes - the default for huge datasets and a metagenome alternative. Unicycler wraps SPAdes, cleans the graph, bridges repeats, and rotates circular replicons to dnaA; hybrid mode (short + ONT/PacBio) is the bacterial finishing standard but a short-only run still cannot beat the repeat limit. ABySS 2.0's Bloom-filter mode enables large-genome assembly on modest hardware but is single-k and niche today.

Pre-Assembly Genome Profiling (run this FIRST)

kmc -k21 -t16 -m64 -ci1 -cs10000 @reads.lst kmcdb tmp/         # k=21: long enough to be mostly unique, short enough for k-mer coverage
kmc_tools transform kmcdb histogram kmer.hist -cx10000
genomescope2 -i kmer.hist -o gs_out -k 21 -p 2                 # -p ploidy; estimates size, heterozygosity, repeat content

GenomeScope2 fits a model to the k-mer frequency histogram and returns, reference-free and before assembly, the genome size, heterozygosity rate, repeat content, and ploidy. A single histogram peak ~ haploid/homozygous; two peaks (a het peak at ~half the homozygous-peak coverage) ~ a heterozygous diploid that will fragment and inflate a short-read assembly. This is the ceiling check: it predicts whether short reads can even produce the assembly being asked for, and it gives the expected size to report assembly total against (assembly >> estimate = un-purged haplotigs, not a big genome).

Per-Method Failure Modes

Adding more coverage to fix contiguity

Trigger: re-sequencing deeper because contigs are short. Mechanism: contiguity is capped by read length vs repeat length, not depth. Symptom: N50 unchanged past ~50-100x; cost wasted. Fix: accept the repeat-determined ceiling, or switch to long reads (-> hifi-assembly).

Hand-picking a single large k for highest N50

Trigger: "we used k=127 because it gave the best N50". Mechanism: large k demands high effective k-mer coverage and is error-fragile; climbing k until N50 peaks is an N50-gaming search that rewards chimeric mega-contigs. Symptom: gappy assembly, or a suspiciously contiguous mis-joined one. Fix: let multi-k auto-select; if overriding, give a small->large list, never a point.