long-read-sequencing/medaka-polishing/SKILL.md
Polishes Oxford Nanopore draft assemblies to higher consensus accuracy with medaka, a basecaller-model-specific neural consensus net, produces haploid variant calls (VCF) for microbial, mitochondrial, or viral samples, and generates amplicon/viral consensus sequences. Covers the model-matching footgun that silently degrades output, why Racon-first is obsolete and medaka runs directly on Flye output as a single pass, why HiFi must never be fed to medaka, the v1->v2 subcommand renames, and the precise medaka_variant deprecation. Use when polishing an ONT-only assembly, generating an amplicon/viral consensus, calling a haploid ONT consensus, or deciding whether medaka, dorado polish, or Clair3 is the right tool.
npx skillsauth add GPTomics/bioSkills bio-long-read-sequencing-medaka-polishingInstall this skill globally with one command. Works with Claude Code, Cursor, and Windsurf.
3 of 9 scanners reported clean
Some scanners were skipped, did not run, or reported a non-clean status. Review each row below.
Reference examples tested with: medaka 2.2+, minimap2 2.28+, samtools 1.19+.
Before using code patterns, verify installed versions match. If versions differ:
<tool> --version then <tool> --help to confirm flagsResults depend on inputs that outlive the binary version - record them:
r1041_e82_400bps_sup_v5.2.0. A mismatch silently degrades output. Prefer auto-detection from the BAM; verify with medaka tools list_models...._sup_v5.2.0, variant ..._sup_variant_v5.0.0 at time of writing); confirm with medaka tools list_models.consensus->inference, stitch->sequence, variant->vcf) and moved the backend to PyTorch; v1 tutorials fail.If code throws an error, introspect the installed tool (medaka --help, medaka_consensus --help) and adapt the example to the actual API rather than retrying.
"Polish my Nanopore assembly" -> Run one medaka consensus pass directly on the assembler output, with the model that matches the basecaller - because a mismatched model silently makes the consensus worse.
medaka_consensus -i reads.fq -d draft.fa -o out/ -t 8 (model auto-detected from the basecaller annotation)medaka is an Oxford Nanopore tool. For PacBio (HiFi/CLR) it is the wrong tool entirely - route to genome-assembly/assembly-polishing.
medaka is a basecaller-model-specific neural consensus net trained on one exact stack (pore + motor enzyme + speed + basecaller mode + basecaller version). Three consequences:
medaka tools resolve_model --auto_model consensus reads.bam); treat a stale model name as a reason to re-basecall, not to proceed.| Mode | Input | medaka's role |
|------|-------|---------------|
| Assembly polishing | Flye/Canu draft + ONT reads | raise per-base QV (homopolymer-indel cleanup is the dominant win) |
| Haploid variant calling | ONT reads + reference (microbial, mito, viral) | medaka_variant wrapper -> haploid VCF (apply with bcftools consensus for a FASTA) |
| Amplicon / viral consensus | tiling-amplicon ONT reads | the non-signal consensus arm of ARTIC fieldbioinformatics / EPI2ME wf-artic |
| Scenario | Recommended | Why |
|----------|-------------|-----|
| ONT-only Flye/Canu assembly | medaka_consensus, ONE pass, auto-detected model | model-matched consensus; racon pre-step is obsolete |
| Native bacterial isolate (modified DNA) | medaka_consensus --bacteria | bacterial-methylation model fixes methylation-motif errors |
| ONT small-variant (diploid/germline) calling | -> clair3-variants | medaka diploid calling deprecated in v2 (Clair3 surpassed it) |
| Haploid microbial/mito/viral VCF | medaka_variant (the renamed haploid wrapper) | still supported in v2 |
| Read-level / human polishing | dorado polish | ONT's emerging successor; identical bacterial weights to medaka today |
| PacBio HiFi/CLR | -> genome-assembly/assembly-polishing | medaka has no PacBio models; never ONT-polish HiFi |
| Unsure which basecaller model produced the reads | re-basecall, then auto-detect | a guessed model silently degrades the consensus |
The wrapper runs three steps: align (mini_align, a thin veil over minimap2 -x map-ont), infer (medaka inference, the neural net over the pileup), and stitch (medaka sequence, regions -> consensus FASTA).
# Canonical modern usage - model auto-detected from the basecaller annotation in the reads
medaka_consensus -i reads.fastq -d draft.fa -o medaka_out/ -t 8
# medaka_out/consensus.fasta is the polished assembly
# Native bacterial isolate: use the methylation-aware bacterial model
medaka_consensus -i reads.fastq -d draft.fa -o medaka_out/ -t 8 --bacteria
# Resolve / list models (do this when auto-detection cannot pick)
medaka tools resolve_model --auto_model consensus reads.bam
medaka tools list_models
medaka runs directly on the assembler (Flye) output as a SINGLE pass - do NOT pre-run Racon (contemporary models are trained on raw assembler output; v2 removed the bundled racon wrapper) and do NOT run medaka twice (iteration was racon's role; a second pass risks flipping correct bases).
medaka_variant emits a VCF only (no consensus FASTA); apply it to the reference with bcftools consensus to get a haploid consensus sequence.
# Wrapper form (renamed from medaka_haploid_variant in v2) - haploid samples only
medaka_variant -i reads.fastq -r reference.fa -o variant_out/
# Manual form - note v2 subcommand names and the hdf -> ref -> out argument order
minimap2 -ax map-ont reference.fa reads.fq | samtools sort -o aln.bam && samtools index aln.bam
medaka inference aln.bam probs.hdf --model r1041_e82_400bps_sup_variant_v5.0.0
medaka vcf probs.hdf reference.fa variants.vcf
# Optional: turn the VCF into a haploid consensus FASTA
bgzip variants.vcf && tabix -p vcf variants.vcf.gz
bcftools consensus -f reference.fa variants.vcf.gz > consensus.fasta
Trigger: running medaka with a model that does not match the basecaller chemistry/version. Mechanism: the net corrects toward the wrong error fingerprint. Symptom: lower held-out QV; medaka exits 0 with no warning. Fix: auto-detect from the BAM; if forced to pick, derive from the actual basecaller and confirm in list_models; treat a stale name as a reason to re-basecall.
Trigger: polishing a PacBio assembly with medaka. Mechanism: ONT-only error model, no PacBio support, on already-QV40+ data. Symptom: degraded/homogenized consensus. Fix: do not; route to genome-assembly/assembly-polishing.
Trigger: running Racon before medaka out of habit. Mechanism: contemporary models are trained on raw assembler output; racon-polished input is off the training distribution. Symptom: no gain or mild harm. Fix: run medaka directly on the Flye output; one pass.
Trigger: an assembly missing a small replicon (~80% identical to a chromosomal region). Mechanism: the absent plasmid's reads misalign onto the chromosome, and medaka "corrects" toward that spurious evidence. Symptom: clustered changes that introduce real errors. Fix: make the assembly structurally complete first; inspect medaka's changes for clustering (clustered = mapping artifact, not scattered homopolymer fixes).
Trigger: judging the polish by medaka's change count or by re-mapping the same reads. Mechanism: medaka optimizes agreement with its input pileup. Symptom: "improvement" that is circular. Fix: reference-free Merqury QV before vs after on held-out / different-platform k-mers.
| Threshold | Source | Rationale | |-----------|--------|-----------| | 1 medaka pass | medaka README | a single trained-model pass; iteration was racon's role, extra passes flip correct bases | | model must match basecaller version | medaka model design | mismatch silently degrades; the #1 ONT-polishing error | | inference threads ~2 | medaka inference behavior | the net is GPU-bound and scales poorly past ~2 CPU threads | | HiFi QV40+ already | EBP/HiFi baseline | nothing for an ONT consensus net to gain; only harm | | measure with held-out Merqury QV | Rhie 2020 | the only honest, reference-free before/after instrument |
| Error / symptom | Cause | Solution |
|-----------------|-------|----------|
| medaka consensus not found / wrong args | v1 subcommand renamed | use medaka inference (or the medaka_consensus wrapper) |
| medaka stitch / medaka variant fail | v1 names | medaka sequence / medaka vcf |
| Polished assembly worse than draft | model mismatch | auto-detect the model; re-basecall if the model is stale |
| medaka errors on PacBio reads | no PacBio models | route to genome-assembly/assembly-polishing |
| Clustered, suspicious changes | missing/mis-structured contig in the draft | complete the assembly first; filter to high-identity alignments |
| Looking for diploid SNP calling | deprecated in v2 | use clair3-variants |
tools
End-to-end CLIP-seq pipeline from FASTQ to ENCODE-compliant binding sites, single-nucleotide crosslink maps, annotation, motifs, and (optionally) differential binding. Use when running the full Yeo lab eCLIP / iCLIP / iCLIP2 / iCLIP3 / irCLIP / PAR-CLIP analysis with SMInput control, protocol-specific UMI extraction, ENCODE STAR parameters, CLIPper or Skipper peak calling with stringent log2 FC and -log10 p thresholds, IDR rescue and self-consistency QC, and downstream motif registration with mCross or PEKA.
development
Detect, date, and contextualize whole-genome duplication (WGD / paleopolyploidy) events using wgd v2 (Chen et al 2024), KsRates (Sensalari 2022 substitution-rate-corrected Ks dating), DupGen_finder (Qiao 2019), MAPS (Li 2018 phylogenomic), POInT (Conant 2008 ordered-block), SLEDGe (2024 ML-based), Whale.jl (Bayesian DL+WGD), and synteny-anchored paranome construction. Use when identifying ancient polyploidy from Ks distributions and synteny block analysis, positioning WGD events relative to speciation, distinguishing tandem from segmental from WGD duplications, dating the 2R/3R vertebrate / fish / salmonid WGDs, building paranome and Ks-age mixture models, applying KsRates substitution-rate correction across lineages, or testing alternative biased-fractionation / dosage-balance models post-WGD.
tools
Build whole-genome alignments using Progressive Cactus (Armstrong 2020 reference-free clade-level WGA), Minigraph-Cactus (Hickey 2024 pangenome-aware), LASTZ chain/net (UCSC pipeline), MUMmer4 (Marçais 2018 pairwise), minimap2 -x asm5/10/20 (Li 2018 fast pairwise), AnchorWave (Song 2022 WGD-aware), and Mauve / progressiveMauve (bacterial). Operates the HAL toolkit (Hickey 2013) for downstream extraction including halSynteny, halLiftover, halBranchMutations, and hal2maf. Use when constructing multi-species alignments for comparative-annotation projection (TOGA), synteny detection, conservation analyses (phyloP / PhastCons), or pangenome graph construction; selecting between reference-free (Cactus) and reference-anchored (LASTZ chains/nets) approaches; tuning sensitivity for closely vs distantly related genomes; or producing HAL files for genome-wide downstream tools.
development
Detect syntenic blocks and structural rearrangements between genomes using MCScanX (Wang 2012), JCVI/MCScan (Tang 2008 Python), GENESPACE (Lovell 2022) for orthology-anchored riparian visualization, SyRI for structural variation, AnchorWave for sequence-level synteny, i-ADHoRe 3.0 for highly diverged species, SynNet for synteny networks, and ntSynt for multi-genome macrosynteny. Use when identifying collinear gene blocks across species, distinguishing macrosynteny from microsynteny, detecting inversions/translocations/duplications, anchoring orthology in WGD lineages, producing publication riparian plots, computing synteny block age via Ks (cross-references whole-genome-duplication), or running synteny-aware ortholog inference in polyploids.