bio-genome-assembly-assembly-qc

What this skill does

## Version Compatibility

Reference examples tested with: QUAST 5.2+, BUSCO 5.5+ (and 6.x for odb12 lineages), compleasm 0.2.6+, Merqury 1.3+, meryl 1.4+, minimap2 2.26+, Inspector 1.2+, CRAQ 1.0+, merfin 1.0+, GenomeScope2 2.0+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags
- Python: `pip show <package>` then `help(module.function)` to check signatures

Results depend on inputs that outlive the binary version - record them:
- BUSCO/compleasm depend on the **lineage dataset** and OrthoDB generation. `_odb10` (BUSCO 5) and `_odb12` (BUSCO 6 default) gene sets are not comparable across the version boundary; a 99% on the shallow `eukaryota_odb10` (~255 genes) is a different claim from 99% on a deep clade set (~5,500+).
- Merqury QV/completeness depend on the **k-mer size** (from `best_k.sh <genome_size>`, not hardcoded) and the **read set** used for the k-mer DB (use accurate reads; see the circularity warning below).
- NG50/NGx/auNG depend on the **expected genome-size estimate** (GenomeScope2 / flow cytometry / a congener).

If code throws an error, introspect the installed tool and adapt rather than retrying.

# Assembly QC

**"Is my genome assembly any good?"** -> Measure all three orthogonal axes - contiguity, completeness, correctness - with reference-free methods, because no single number (least of all N50) is quality.
- CLI: `quast.py asm.fa --large --eukaryote -o out` (contiguity + reference-based structure), `busco -i asm.fa -m genome -l <lineage>` or `compleasm run -a asm.fa -l <lineage>` (gene completeness), `merqury.sh reads.meryl asm.fa out` (reference-free QV + k-mer completeness), `inspector.py -c asm.fa -r reads.fq` (reference-free structural errors)

## The Single Most Important Modern Insight -- Quality Is Three Orthogonal Axes; N50 Is the Most-Gamed One

Assembly quality is **three genuinely orthogonal axes - contiguity, completeness, correctness - and a single number on any one is not quality.** The axes do not predict each other, and the diagnostic failure modes prove it:

- **Contiguous + wrong:** a single-contig "chromosome" that is three chromosomes misjoined. Perfect N50, catastrophic correctness. Only Hi-C / a same-species reference / read-discordance catches it.
- **Complete + shredded:** BUSCO 99%, but repeats collapsed, segmental duplications merged, intergenic space wrong. BUSCO is gene-space-only and cannot see it.
- **Accurate + incomplete:** QV60 over the 92% that assembled, with the hard 8% (centromeres, rDNA, satellites) simply absent. QV is silent about what is not there.

The field's historical sin is reporting **contiguity alone** because it is cheapest to compute and easiest to game. **N50 is the most-gamed metric in genomics:** it rises when sequence is thrown away (N50 is computed on what survives), when misjoins are *not* broken (a misjoined contig is a long contig), and when haplotigs are retained. A bigger N50 is louder, not better. Three load-bearing moves:

1. **Report auN/NGx, not bare N50.** auN = the area under the Nx curve = length-weighted mean contig length; it integrates the whole curve and is continuous where N50 jumps discontinuously (the small-L50 / T2T regime). NGx/auNG normalize to the **expected genome size**, coupling contiguity to completeness (an assembly that drops half the genome gets a great N50 but a terrible NG50). Always report contig AND scaffold N50 - if scaffold >> contig, the contiguity is glue (Ns), not sequence.
2. **Default to reference-free.** For a *novel* genome there is no trusted reference; QUAST against a divergent relative reports real inversions/SVs as "misassemblies" and real SNPs as "mismatches". Use Merqury QV (accuracy) + Inspector/CRAQ (structure) + asmgene (false dup/collapse). QUAST is the special case "I have a same-organism reference," not the default.
3. **Report a Merqury QV - and never compute it on the polishing reads.** QV is the reference-free accuracy standard reviewers now demand; an assembly paper with no QV is a red flag. But QV from the same reads used for polishing is **circular** - the polisher already made the assembly agree with those reads, so the QV measures convergence, not correctness. Build the k-mer DB from accurate, ideally independent reads (HiFi/Illumina, not noisy ONT).

## Tool Taxonomy

| Tool | Citation | Axis / Role | When |
|------|----------|-------------|------|
| QUAST / calN50 | Gurevich 2013 *Bioinformatics*; auN = Li blog (no journal) | contiguity (auN/NGx, N50/L50) + reference-based structure (NA50, misassemblies) | always for contiguity; structure only vs a same-organism reference |
| BUSCO | Manni 2021 *Mol Biol Evol*; Simão 2015 *Bioinformatics* | gene-space completeness (C/S/D/F/M) | universal; conservative on good genomes |
| compleasm | Huang & Li 2023 *Bioinformatics* | gene-space completeness, miniprot-based | faster + more sensitive; default on HiFi/T2T-era genomes |
| Merqury | Rhie 2020 *Genome Biol* | reference-free QV + k-mer completeness + spectra-cn + phasing | always; the accuracy gold standard |
| merfin | Formenti 2022 *Nat Methods* | multiplicity-corrected QV / polishing | refine QV biased by k-mer multiplicity |
| Inspector | Chen 2021 *Genome Biol* | reference-free structural + base errors (long reads) | novel genomes; can also correct |
| CRAQ | Li 2023 *Nat Commun* | reference-free structural/regional errors (clipped alignments) | novel genomes; flags misjoins to split |
| asmgene | Li (minimap2, no separate journal) | gene collapse / false duplication | high-quality genomes where BUSCO saturates |
| GenomeScope2 | Ranallo-Benavidez 2020 *Nat Commun* | genome size / het / repeat % from k-mers | the size estimate NG50/auNG/spectra-cn need (-> genome-profiling) |

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Novel genome, no trusted reference | Merqury QV + k-mer completeness + BUSCO/compleasm + auN/NGx + Inspector/CRAQ | reference-free triad; QUAST structure is uninterpretable here |
| Same-species (near-isogenic) reference available | add QUAST `--large --eukaryote -r ref.fa` (NA50, misassemblies) | reference-based structure is trustworthy only vs the same organism |
| High-quality HiFi/T2T-era assembly, BUSCO looks low | compleasm | BUSCO under-reports good genomes (its predictor, not the assembly, misses genes) |
| Contiguity claim must resist gaming | auN/auNG via `calN50.js -L <size>` | N50 is a single unstable order-statistic and is gameable |
| High BUSCO-Duplicated, size > expected | spectra-cn + asmgene + GenomeScope2 size -> purge_dups | distinguish uncollapsed haplotigs (purge) from real WGD (keep) |
| Phased diploid / trio assembly | Merqury hap-mers: switch/hamming error + blob plot | phasing accuracy is the extra axis |
| Need genome size / het before NG50 | -> genome-profiling (GenomeScope2) | NG/auN/spectra-cn all need a size estimate |
| Reads not yet QC'd | -> read-qc/quality-reports | garbage-in caps assembly quality |
| Bacterial isolate / MAG completeness+contamination | -> contamination-detection (CheckM2/GUNC/MIMAG) | marker-gene completeness/contamination is a different problem |

## Contiguity -- auN/NGx (not bare N50)

```bash
k8 calN50.js -L <genome_size> asm.fa       # N50/L50 + NG50/NGx + auN/auNG; -L sets genome size for NG/auNG (ships with minimap2)
quast.py asm.fa --large --eukaryote -t 16 -o quast_out   # N50/L50, GC, # contigs; NG50 only with -r or --est-ref-size; structure only if -r given
```

`--large` implies `--eukaryote --min-contig 3000 --min-alignment 500 --extensive-mis-size 7000`. Report contig AND scaffold N50; a scaffold N50 far above the contig N50 means the contiguity is scaffolding Ns, and every gap is a join hypothesis that could be a misassembly. NA50 (QUAST, contigs broken at misassemblies) far below N50 means the contiguity is partly fictional.

## Completeness -- BUSCO / compleasm + Merqury k-mer c