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)development
Installs 425 bioinformatics skills covering sequence analysis, RNA-seq, single-cell, variant calling, metagenomics, structural biology, and 56 more categories. Use when setting up bioinformatics capabilities or when a bioinformatics task requires specialized skills not yet installed.
testing
Chains a somatic (tumor-normal) SNV/indel and structural-variant pipeline end to end with GATK Mutect2 (or Strelka2), wiring the somatic-specific machinery - panel-of-normals and gnomAD germline-resource priors, GetPileupSummaries/CalculateContamination, and LearnReadOrientationModel FFPE/oxoG orientation-bias filtering fed into FilterMutectCalls. Use when calling somatic mutations from a tumor-normal pair (or tumor-only with PoN caveats), deciding which artifact filter removes which class of false positive, reasoning about VAF/purity/ploidy and clonal-vs-subclonal detection, adding somatic SV/CNV or TMB/MSI/signatures, or routing variants to AMP/ASCO/CAP tier and oncogenicity interpretation (never germline ACMG).
development
End-to-end pooled and single-cell CRISPR screen analysis from FASTQ to hit genes. Orchestrates library design QC, guide counting, six-stage screen QC (plasmid Gini, replicate Pearson, CEGv2 PR-AUC, copy-number artifact), method-appropriate hit calling across MAGeCK RRA/MLE, BAGEL2, drugZ, JACKS, and Chronos, cancer-cell-line copy-number correction (CRISPRcleanR / Chronos), batch correction for multi-batch screens, and the specialized branches for combinatorial paralog screens, single-cell Perturb-seq, base-editor variant-function screens, prime-editor screens, and in vivo bottleneck-aware screens. Use when analyzing any pooled CRISPR screen end-to-end, matching the hit-calling method to the experimental design, integrating copy-number correction into the pipeline, or branching the workflow for single-cell, combinatorial, base-editor, prime-editor, or in vivo variants.
development
Transcribe DNA to RNA and translate to protein using Biopython, with NCBI codon-table selection, CDS validation, and six-frame ORF finding. Use when converting a CDS or ORF to its amino-acid sequence, selecting a non-standard (mitochondrial, bacterial, ciliate) genetic code, validating a coding sequence, or scanning all reading frames.