bio-genome-assembly-assembly-polishing

What this skill does

## Version Compatibility

Reference examples tested with: Racon 1.5+, medaka 2.0+, minimap2 2.26+, bwa 0.7.17+, samtools 1.19+, Polypolish 0.6+, pypolca 0.3+, Pilon 1.24+, Merqury 1.3+, meryl 1.4+.

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

The **medaka consensus model** matters more than the medaka binary version: the model must match the basecaller + pore chemistry + caller mode + version (`-m r1041_e82_400bps_sup_v5.0.0`); a mismatched model silently degrades the consensus. List options with `medaka tools list_models`. medaka's status is a moving target - ONT now steers toward `dorado polish` and has deprecated medaka's diploid-variant workflow in favour of Clair3; verify current guidance and the `medaka_consensus` wrapper vs newer subcommands. If code throws an error, introspect the installed tool and adapt rather than retrying.

# Assembly Polishing

**"Polish my genome assembly"** -> Decide whether base-level consensus correction is warranted, then apply a read-type-matched polisher and measure the gain with held-out k-mers - or, for an already-accurate assembly, decline.
- CLI: `racon reads.fq aln.sam draft.fa` then `medaka_consensus -i reads.fq -d draft.fa -o out -m <model>` (ONT); `polypolish polish draft.fa f1.sam f2.sam` + `pypolca run` (hybrid/short); `merqury.sh reads.meryl asm.fa out` (measure QV before/after)

## The Single Most Important Modern Insight -- Polishing Is Being Engineered Out, and Reflexive Polishing Does Net Harm

Polishing exists to fix the characteristic error of *noisy* long reads: indels in homopolymers and low-complexity tracts (the pore/RT stutters on `AAAAAA`) plus residual substitutions. As read accuracy rose (PacBio HiFi ~Q40-50 reads, ONT R10.4.1 + Dorado `sup`/duplex), the assembly consensus is already near-perfect, and a mapping-based polisher's read-pileup step injects more errors than it removes. Two load-bearing rules dominate:

1. **Do NOT polish HiFi assemblies by default.** A hifiasm/HiCanu HiFi assembly starts at ~Q40+; mapping-based polishers (Racon, Pilon, GCpp) routinely *lower* QV and introduce haplotype-switch errors. Only polish HiFi with a tool explicitly built not to overcorrect (NextPolish2, DeepPolisher) and only if Merqury QV says there is a real deficit. For HiFi the burden of proof is on polishing, not against it.
2. **Measure the gain with reference-free, HELD-OUT Merqury k-mers - never the reads polished with, never BUSCO.** A polisher's literal objective is to maximize agreement with its input reads, so scoring it with those same reads is circular and always looks like an improvement. The honest stop signal is a Merqury QV *plateau*, not a fixed iteration count - and polishing can drive QV *down* silently. Polish with one platform, evaluate with another (polish with ONT, measure with Illumina k-mers).

Polishing is transitional infrastructure: indispensable in the CLR/early-ONT error era, shrinking toward irrelevance as raw accuracy climbs upstream (better basecalling, hifiasm's own consensus). The skill's job is to say *when not to polish* as firmly as how to.

## Tool Taxonomy

| Tool | Citation | Role | When |
|------|----------|------|------|
| Racon | Vaser 2017 *Genome Res* | POA read-to-assembly consensus; fast workhorse | first-pass ONT/CLR self-polish (1-4 rounds); needs minimap2 SAM/PAF |
| medaka | nanoporetech/medaka (software) | ONT neural-network consensus | one pass *after* Racon; model MUST match basecaller |
| dorado polish | ONT/dorado (software) | modern ONT consensus, owns basecalling context | the emerging medaka replacement; verify current status |
| Polypolish | Wick & Holt 2022 *PLoS Comput Biol* | Illumina polish using ALL alignments (`bwa mem -a`) | repeat-rich long-read assemblies; best-in-class short-read polish |
| pypolca | Bouras 2024 *Microb Genom* (POLCA: Zimin 2020) | fast Illumina SNP/indel correction | pair with Polypolish (modern bacterial practice) |
| Pilon | Walker 2014 *PLoS ONE* | all-in-one SNP/indel/local short-read fix | legacy; small genomes only; mismaps in repeats, ~1 GB heap/Mb |
| NextPolish / NextPolish2 | Hu 2020 *Bioinformatics* / Hu 2024 *GPB* | iterative LR+SR polish / HiFi repeat+phase-aware | NextPolish2 is the HiFi-safe successor |
| Hapo-G | Aury & Istace 2021 *NARGAB* | haplotype-aware short-read polish | heterozygous diploids; preserves both alleles |
| ntEdit | Warren 2019 *Bioinformatics* | Bloom-filter k-mer polish, no mapping | very large (>3 Gb) genomes; scales where alignment can't |
| DeepPolisher | Mastoras 2025 *Genome Res* | transformer correction (PHARAOH ONT phasing) | HiFi/T2T human; halves errors without overcorrection |
| PEPPER-Margin-DeepVariant | Shafin 2021 *Nat Methods* | haplotype-aware ONT consensus via variant calling | human ONT-only polishing |
| DeepConsensus | Baid 2023 *Nat Biotechnol* | improves CCS *reads* UPSTREAM of assembly | NOT a polisher - operates before assembly (common category error) |

## Decision Tree by Scenario

| Scenario (data type) | Recommended | Why |
|----------------------|-------------|-----|
| PacBio HiFi only | **Do NOT polish by default** -> check Merqury QV first | already ~Q40-50; Racon/Pilon lower QV + collapse haplotypes |
| HiFi with a real Merqury QV deficit | NextPolish2 or DeepPolisher | repeat/phase-aware; won't overcorrect het sites |
| ONT-only, noisy (R9 / R10 hac) | Racon (1-4 rounds) -> medaka (1 pass), or dorado polish | canonical ONT chain; neural consensus on the ONT error model |
| ONT-only, human/diploid | PEPPER-Margin-DeepVariant or DeepPolisher | haplotype-aware; preserves both alleles |
| Long-read draft + Illumina (hybrid) | long-read pass, then Polypolish + pypolca | SR fixes residual homopolymer indels Racon/medaka missed |
| Short-read-only assembly (SPAdes) | usually no polishing needed | Illumina is already ~Q40+; structural sins survive anyway |
| HiFi + Illumina | use Illumina k-mers for *evaluation*, not correction | Merqury hybrid QV, not a polishing pass |
| Heterozygous / repeat-rich genome | Polypolish (`-a`), Hapo-G, NextPolish2 - or none | non-haplotype-aware polishers erase het / homogenize paralogs |
| Reads not yet QC'd / wrong basecaller | -> read-qc/quality-reports, -> long-read-sequencing/long-read-alignment | garbage-in or model-mismatch makes polishing worse than skipping |
| Complaint is "fragmented" not "low QV" | -> not polishing (scaffolding/gap-filling) | polishing fixes bases, not contiguity |
| Map reads before polishing | short -> read-alignment/bwa-alignment; long -> long-read-sequencing/long-read-alignment | polishers consume a BAM/SAM/PAF, not raw reads |

## The Canonical ONT Chain: Racon -> medaka

**Racon** does the bulk cheap consensus from the long reads themselves; it is a standalone consensus module and does NOT map - the SAM/PAF is supplied. Re-map every round (1-4 rounds; gains decay fast).

```bash
minimap2 -t 16 -ax map-ont draft.fasta reads.fq.gz > aln.sam   # map-pb for PacBio CLR
racon -t 16 reads.fq.gz aln.sam draft.fasta > racon1.fasta
# re-map reads to racon1.fasta and repeat for round 2... measure QV each round, stop at plateau
```

**medaka** applies an ONT-trained neural model for the final consensus - exactly ONE pass after the Racon rounds (it is not an iterate-many-times tool; running medaka twice is a tell). For medaka mechanics and model tables see long-read-sequencing/medaka-polishing; this skill owns the *strategy*.

```bash
medaka_consensus -i reads.fq.gz -d racon_final.fasta -o medaka_out -t 16 \
  -m r1041_e82_400bps_sup_v5.0.0          # model MUST match basecaller+chemistry+caller+version
```

Recent basecallers embed the model in the FASTQ so medaka auto-selects; if the data was basecalled with an old/unknown caller, pick manually from `medaka tools list_models`, and treat a stale/deprecated model name as a reason to rebasecall rather than proceed.

## Hybrid / Short-Read Polishi