Bio Genome Assembly Assembly Polishing

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

Decision Tree by Scenario

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).

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.

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 Polishing: Polypolish + pypolca

Polypolish aligns short reads to all locations (bwa mem -a) so it can disambiguate which repeat copy a read belongs to instead of forcing one placement - this is how it beats Pilon in repeats and why it almost never introduces errors.

bwa index draft.fasta
bwa mem -t 16 -a draft.fasta reads_1.fq.gz > aln_1.sam     # -a = ALL alignments (the whole point)
bwa mem -t 16 -a draft.fasta reads_2.fq.gz > aln_2.sam
polypolish filter --in1 aln_1.sam --in2 aln_2.sam --out1 filt_1.sam --out2 filt_2.sam
polypolish polish draft.fasta filt_1.sam filt_2.sam > polypolish.fasta   # --careful (v0.6+) for low depth
pypolca run -a polypolish.fasta -1 reads_1.fq.gz -2 reads_2.fq.gz -o pypolca_out -t 16 --careful

Bouras 2024 depth-tiered bacterial recommendation: depth <5x -> Polypolish --careful alone; 5-25x -> Polypolish --careful + pypolca --careful; >25x -> Polypolish (default) + pypolca --careful. Modern bacterial best practice is Polypolish + pypolca, NOT Pilon.

Pilon (Legacy Short-Read)

bwa mem -t 16 draft.fasta r1.fq r2.fq | samtools sort -o frags.bam
samtools index frags.bam
java -Xmx16G -jar pilon.jar --genome draft.fasta --frags frags.bam --output pilon --fix all --changes

--fix modes: snps, indels, bases (=snps+indels), gaps, local, all (default), none. Pilon needs a sorted+indexed BAM (route mapping to read-alignment/bwa-alignment). It is legacy: it uses best-placement alignments so it mismaps in repeats (miscorrects toward paralogs), and its ~1 GB heap per Mb of genome OOMs on large eukaryotes - both reasons Polypolish/pypolca displaced it.

Measuring the Gain: Held-Out Merqury QV

Goal: Decide whether a polish actually helped, without fooling yourself.

Approach: Build a meryl k-mer DB from an independent / different-platform read set (k from Merqury's best_k.sh, not hardcoded), then run Merqury on the pre- and post-polish assemblies and compare QV. A polish that does not raise QV did not help; one that lowers it must be reverted.

K=$(sh $MERQURY/best_k.sh 5000000 | tail -n1 | awk '{print int($1+0.5)}')   # K from genome size; round float->int
meryl count k=$K output reads.meryl illumina_reads.fq.gz               # eval reads != polishing reads
merqury.sh reads.meryl draft.fasta    qv_before                        # QV of the input
merqury.sh reads.meryl polished.fasta qv_after                         # QV after polishing

Use a held-out set or a different platform than was polished with (polish with ONT, evaluate with Illumina k-mers); the CHM13 effort triangulated with both HiFi and Illumina k-mers (Mc Cartney 2022). Do NOT use BUSCO as the polishing metric - it measures gene-space completeness and barely moves with the homopolymer-indel QV that polishing changes, so a flat BUSCO masks a silent QV drop. A per-gene internal-stop / frameshift count is a useful secondary readout. See assembly-qc for the full QV/spectra-cn workflow.

Per-Method Failure Modes