long-read-sequencing/clair3-variants/SKILL.md
Calls germline small variants (SNPs and indels) from Oxford Nanopore and PacBio HiFi long reads with Clair3, a two-stage (pileup + full-alignment) deep-learning caller, selecting the chemistry- and basecaller-version-matched model, enabling read-based phasing, and benchmarking against GIAB with stratification. Covers why the model string is the experiment (no auto-detection, silent degradation on mismatch), why ONT homopolymer/STR indels are the residual error whole-genome F1 hides, and the somatic/trio/RNA boundary to the ClairS/Clair3-Trio family. Use when calling germline SNVs/indels from ONT or HiFi BAMs, choosing a Clair3 model, phasing variants, or benchmarking long-read calls.
npx skillsauth add GPTomics/bioSkills bio-long-read-sequencing-clair3-variantsInstall 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: Clair3 2.0+, whatshap 2.0+, bcftools 1.19+, hap.py 0.3.15+.
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_v500). There is NO auto-detection; --model_path is mandatory and a mismatch silently degrades calls.pileup.pt/full_alignment.pt. v1 TensorFlow models do NOT load in v2._with_mv signal-aware) lives in the rerio clair3_models/ repo; only a subset is bundled.If code throws an error, introspect the installed tool (run_clair3.sh --help) and adapt the example to the actual API rather than retrying.
"Call variants from my long reads" -> Run Clair3 with the model that matches how the reads were basecalled, phase, and benchmark with stratification - because the model string, not the command, determines accuracy.
run_clair3.sh --bam_fn=aln.bam --ref_fn=ref.fa --output=out/ --threads=16 --platform=ont --model_path=/models/r1041_e82_400bps_sup_v500Scope: germline diploid SNPs + small indels. NOT structural variants (-> structural-variants), NOT somatic/mosaic (-> ClairS/ClairS-TO), NOT RNA (-> Clair3-RNA).
Clair3's accuracy is gated by two facts a naive user misses:
--model_path at a specific model folder. Three axes must ALL match: chemistry (r941 vs r1041), basecaller tier (fast/hac/sup), and basecaller version (g5014/v430/v500/v520), plus the optional _with_mv signal-aware axis if the BAM has Dorado mv tags. Wrong model = no crash, no warning, measurably worse calls (indels most). Derive the model from the basecaller string in the run metadata; pick the model version closest to but not above the basecaller version.Clair3 "symphonizes" two networks: a fast pileup model (summarized per-position statistics) that calls the large majority of sites, and a slow full-alignment model (haplotype-resolved read tensor) that re-evaluates only the uncertain subset. Internally Clair3 phases the top het-SNP pileup calls with WhatsHap, haplotags the BAM, and feeds the haplotagged reads to the full-alignment model - which is why read-based phasing buys ~6% indel F1, not cosmetics. Output: merge_output.vcf.gz (final).
Model name anatomy (r1041_e82_400bps_sup_v500): pore (r1041=R10.4.1), flowcell (e82), speed (400bps), basecaller tier (sup/hac/fast), basecaller version (v500=Dorado 5.0.0, g5014=Guppy 5.0.14). The _with_mv suffix uses Dorado move-table tags for best accuracy when present.
| Data | --platform | Model |
|------|--------------|-------|
| ONT R10.4.1 sup, Dorado v5.x, mv tags present | ont | r1041_e82_400bps_sup_v520_with_mv |
| ONT R10.4.1 sup, Dorado v5.0.0 | ont | r1041_e82_400bps_sup_v500 |
| ONT R10.4.1 hac | ont | r1041_e82_400bps_hac_v500/_v520 |
| ONT R9.4.1 (any tier) | ont | r941_prom_sup_g5014 |
| PacBio HiFi Revio | hifi | hifi_revio |
| PacBio HiFi Sequel II | hifi | hifi_sequel2 |
| Illumina (supported) | ilmn | ilmn |
| PacBio CLR | - | not supported -> PEPPER-Margin-DeepVariant |
| Scenario | Tool | Why |
|----------|------|-----|
| Germline SNV/indel, single sample | Clair3 | this skill |
| Somatic, paired tumor-normal | ClairS | VAF-aware; Clair3 germline priors cannot find low-VAF somatic |
| Somatic, tumor-only | ClairS-TO | tumor-only ensemble |
| De novo / Mendelian trio | Clair3-Nova / Clair3-Trio | family-aware |
| Long-read RNA variants | Clair3-RNA | RNA model |
| ONT R10.4.1, also considering DeepVariant | either | neck-and-neck on R10 sup; native-ONT DeepVariant (Kolesnikov 2024) superseded PEPPER-Margin |
| Non-human / draft / bacterial reference | Clair3 + --include_all_ctgs | default calls only chr1-22,X,Y -> empty output otherwise |
| Cohort joint genotyping | Clair3 gVCF -> GLnexus | bcftools merge on gVCFs is NOT joint genotyping |
# Germline ONT calling (model MUST match the basecaller)
run_clair3.sh \
--bam_fn=aln.bam --ref_fn=ref.fa --output=clair3_out/ \
--threads=16 --platform=ont \
--model_path=/opt/models/r1041_e82_400bps_sup_v500
# final VCF: clair3_out/merge_output.vcf.gz
# Phase the final output VCF (WhatsHap); --longphase_for_phasing swaps only the INTERNAL
# phaser to LongPhase (faster, SV-aware). For a LongPhase-phased final VCF use
# --use_longphase_for_final_output_phasing instead of --enable_phasing.
run_clair3.sh ... --enable_phasing --longphase_for_phasing
# Phased calls go to clair3_out/phased_merge_output.vcf.gz; merge_output.vcf.gz stays UNPHASED.
# Non-human / draft assembly reference - call ALL contigs
run_clair3.sh ... --include_all_ctgs
# Targeted / amplicon panel
run_clair3.sh ... --bed_fn=panel.bed --gvcf
# Benchmark against GIAB with stratification (the step that reveals ONT indel errors)
hap.py giab_truth.vcf.gz clair3_out/merge_output.vcf.gz \
-f giab_confident.bed -r ref.fa --engine=vcfeval \
--stratification giab_stratifications.tsv -o bench/hg002
Trigger: --model_path pointing at a model that does not match the basecaller chemistry/tier/version. Mechanism: no auto-detection; the wrong network runs. Symptom: no error, lower F1 (indels most). Fix: derive the model from the basecaller string; verify the folder exists (rerio for the full set); for v2 ensure .pt models.
Trigger: bacterial genome or draft assembly without chr1-22,X,Y names. Mechanism: Clair3 calls only standard human contigs by default. Symptom: near-empty VCF. Fix: --include_all_ctgs.
Trigger: reporting only whole-genome F1. Mechanism: ONT indel errors concentrate in homopolymer/STR/low-complexity strata. Symptom: ~99.5% global indel F1 but much lower in LowComplexity. Fix: stratify with GIAB BEDs (Dwarshuis 2024); use CMRG for medically relevant genes.
Trigger: lowering --snp_min_af/--indel_min_af to catch low-VAF variants. Mechanism: germline model expects ~0.5/1.0 allele fractions, is not VAF-aware. Symptom: germline-model false positives at low AF, missed true somatic. Fix: ClairS (paired) / ClairS-TO (tumor-only).
Trigger: an old TensorFlow model dir with Clair3 v2. Mechanism: v2 needs PyTorch .pt models. Symptom: model load failure. Fix: use pileup.pt/full_alignment.pt models (Converted Rerio).
| Threshold | Source | Rationale |
|-----------|--------|-----------|
| Recommended depth ~20-60x | Clair3 guidance | sensitivity (hets, indels) falls off below ~20x; --min_coverage default 2 is a floor, not a recommendation |
| Phasing buys ~6% indel F1 | Zheng 2022 | haplotagged reads disambiguate indel alleles in repeats |
| ONT R10.4.1 sup: SNP F1 ~99.99%, indel F1 ~99.5% | GIAB benchmarks | indel residual lives in homopolymer/STR strata |
| --var_pct_full 0.3 (default) | Clair3 README | fraction of low-quality pileup calls re-run by full-alignment; raise for recall, slower |
| Stratify with GIAB / CMRG | Dwarshuis 2024 | global F1 hides the ONT indel problem |
| Error / symptom | Cause | Solution |
|-----------------|-------|----------|
| Empty/near-empty VCF on non-human ref | default calls only chr1-22,X,Y | --include_all_ctgs |
| Model fails to load | v1 TF model with v2 Clair3 | use .pt (PyTorch) models |
| --model_path .../models/ont not found | no generic ont/hifi model | point at a specific model subfolder |
| Worse-than-expected indels | wrong-version or wrong-tier model | match the basecaller model exactly |
| "joint genotyping" gave odd merges | bcftools merge on gVCFs is not joint calling | use GLnexus |
| Looking for somatic/low-VAF variants | germline caller | use ClairS / ClairS-TO |
--MD; use minimap2 >=2.28)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.