name: bio-variant-calling description: Call germline SNPs and indels from a BAM/CRAM with bcftools mpileup and call, and select the right calling engine for the job. Use when generating a VCF from aligned reads, choosing between bcftools, GATK HaplotypeCaller, DeepVariant, and DRAGEN, setting ploidy for haploid/organelle/polyploid/sex-chromosome calling, or deciding whether pileup-based calling is good enough versus a local-reassembly caller for indels and difficult regions. Not for cohort joint genotyping (see variant-calling/joint-calling), GATK-specific workflows (see variant-calling/gatk-variant-calling), deep-learning calling (see variant-calling/deepvariant), or somatic/low-VAF detection. origin: openai4s category: bioskills/variant-calling metadata: tool_type: cli primary_tool: bcftools third_party: name: GPTomics/bioSkills repository: https://github.com/GPTomics/bioSkills commit: d91ed3d563019e649dc854c56ccd62551359488a license: MIT
Reference examples tested with: bcftools 1.19+
Before using code patterns, verify installed versions match. If versions differ:
- CLI: <tool> --version then <tool> --help to confirm flags
If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
Note: bcftools mpileup applies BAQ (per-Base Alignment Quality) by default; this is a real behavior that changes calls, not a nuisance flag (see The Governing Principle).
"Call SNPs and indels from my aligned reads" -> Compute per-position genotype likelihoods from a BAM/CRAM against the reference, then call variant sites under a Bayesian model at the assumed ploidy.
- CLI (fast, position-based): bcftools mpileup -f ref.fa in.bam | bcftools call -mv
- CLI (reassembly, higher indel accuracy): GATK HaplotypeCaller (variant-calling/gatk-variant-calling)
- CLI (deep learning): DeepVariant (variant-calling/deepvariant)
This skill does the bcftools calling and is the engine-selection hub: it tells the agent when pileup calling is the right tool and when to hand off to a reassembly or deep-learning caller.
There are two families of short-variant caller, and the choice between them is the single most consequential decision here.
Consequence: bcftools is fine-to-excellent for simple germline SNPs and quick genome-wide scans, materially weaker on indels and in low-complexity / segmental-duplication / MHC regions, and not built for somatic low-VAF detection or scalable cohort joint calling. Pick the engine from the analysis, not from habit.
Guidance, not dogma; on a production human pipeline, validate against current GIAB/GA4GH benchmarks (hap.py + vcfeval) before committing.
| Engine | Best when | Fails / weak when | Hand off to |
|---|---|---|---|
| bcftools mpileup|call | Simple germline SNPs; non-model/organelle/microbial genomes (no training data, any ploidy); quick exploratory scans; low compute; small multi-sample sets | Indels in homopolymers/STRs; segdups, MHC, low-mappability; low-VAF somatic/mosaic; cohorts beyond ~100 samples | this skill |
| GATK HaplotypeCaller | Auditable open-source human WGS/WES; every parameter inspectable; the joint-calling/best-practices orthodoxy (GVCF -> GenomicsDBImport -> GenotypeGVCFs) | Lower indel/difficult-region accuracy than DeepVariant/DRAGEN; local assembly can abort in pathological high-depth/repeat regions | variant-calling/gatk-variant-calling; variant-calling/joint-calling |
| DeepVariant | Best open-source accuracy on indels and difficult regions; PacBio HiFi / ONT (platform-specific trained models); generalizes off one training sample | Needs the correct platform model (wrong model degrades accuracy); GPU helps; cohort merge needs GLnexus, not GenotypeGVCFs | variant-calling/deepvariant |
| DRAGEN | Maximum throughput on Illumina (FPGA, ~20-25 min/genome); leads difficult-to-map benchmarks (alt-aware mapping) | Proprietary/hardware- or license-gated; ML recalibrator trained on GIAB truth (benchmark-overfitting caveat) | vendor pipeline; HaplotypeCaller --dragen-mode for an open-source approximation |
Honest state of the field: DeepVariant and DRAGEN lead on indels and difficult regions; GATK is the joint-calling and best-practices reference everyone else is measured against; bcftools wins on speed, simplicity, non-model organisms, and organelle/haploid calling. On easy SNPs every modern caller exceeds F1 0.999, so a caller's headline SNP number is rarely the deciding factor - indels and hard regions are.
Goal: Detect germline SNPs and indels from aligned reads with the pileup-and-call pipeline.
Approach: Generate per-position genotype likelihoods with mpileup (BAQ on by default), pipe as uncompressed BCF into the multiallelic caller.
bcftools mpileup -f reference.fa input.bam | bcftools call -mv -Oz -o variants.vcf.gz
bcftools index variants.vcf.gz
# -Ou between steps avoids VCF (de)serialization; -q/-Q drop poorly-supported reads/bases;
# -a requests the FORMAT tags downstream filtering needs (DP, allelic depths, strand-bias p)
bcftools mpileup -Ou -f reference.fa \
-q 20 -Q 20 \
-a FORMAT/DP,FORMAT/AD,FORMAT/SP \
input.bam | \
bcftools call -mv -Oz -o variants.vcf.gz
bcftools index variants.vcf.gz
# Single region / BED targets
bcftools mpileup -f reference.fa -r chr1:1000000-2000000 input.bam | bcftools call -mv -Oz -o region.vcf.gz
bcftools mpileup -f reference.fa -R targets.bed input.bam | bcftools call -mv -Oz -o targets.vcf.gz
# Multiple BAMs (small cohorts only; see The Governing Principle for the scaling limit)
bcftools mpileup -f reference.fa sample1.bam sample2.bam sample3.bam | bcftools call -mv -Oz -o cohort.vcf.gz
# BAM list file: one path per line
bcftools mpileup -f reference.fa -b bams.txt | bcftools call -mv -Oz -o cohort.vcf.gz
| Stage | Flag | Effect |
|---|---|---|
| mpileup | -f ref.fa |
Reference FASTA (required); must be the exact one used for alignment |
| mpileup | -q INT |
Min mapping quality; -q 20 drops ambiguously placed reads (paralog mismapping) |
| mpileup | -Q INT |
Min base quality; -Q 20 drops low-confidence base calls |
| mpileup | -a LIST |
Extra FORMAT/INFO tags: FORMAT/AD (allelic depths), FORMAT/DP, FORMAT/SP (Phred strand-bias p), FORMAT/ADF/ADR (per-strand), INFO/AD |
| mpileup | -d INT |
Max per-file depth (default 250); set to 3-4x expected mean coverage to avoid truncating high-coverage sites |
| mpileup | -B / -E |
-B disables BAQ (more raw indel signal, more false SNPs near indels); -E recomputes BAQ on the fly (more sensitive, slower) |
| call | -m |
Multiallelic caller - default, recommended for all new work |
| call | -c |
Consensus caller - legacy; only for reproducing old pipelines |
| call | -v |
Emit variant sites only (omit to emit all sites, e.g. for hom-ref confidence) |
| call | -O z\|b\|u\|v |
Output: z bgzipped VCF, b BCF, u uncompressed BCF (piping), v VCF |
| call | --ploidy / --ploidy-file |
Sample/region ploidy (below) |
| call | -P FLOAT |
Mutation-rate prior (default 1.1e-3, human); lower for inbred lines, raise for diverse/outbred populations |
The multiallelic caller (-m) handles sites with several ALT alleles natively and is statistically superior; the consensus caller (-c) exists only for backward reproducibility.
Goal: Match the caller's ploidy to the biology so genotypes are representable.
Approach: Set a scalar ploidy for uniform samples, or a ploidy file (or built-in preset) to vary ploidy by region and sex.
Wrong ploidy silently corrupts calls: calling a diploid as haploid halves heterozygous sensitivity; calling a haploid/hemizygous region as diploid manufactures false heterozygous calls from every error and paralog mismap.
# Haploid: bacteria, mitochondria (nuclear germline heteroplasmy caveat below), non-PAR chrX/chrY in a male
bcftools mpileup -f reference.fa input.bam | bcftools call -m --ploidy 1 -Oz -o haploid.vcf.gz
# Built-in human preset applies karyotype-aware sex-chromosome ploidy
bcftools call -m --ploidy GRCh38 ...
# Ploidy file: CHROM FROM TO SEX PLOIDY (chrY absent in females -> 0)
# chrX 1 -1 M 1
# chrX 1 -1 F 2
# chrY 1 -1 M 1
# chrY 1 -1 F 0
# * 1 -1 * 2
bcftools mpileup -f reference.fa input.bam | bcftools call -m --ploidy-file ploidy.txt -Oz -o sexaware.vcf.gz
Scope notes: true mitochondrial heteroplasmy is continuous-VAF (not 0/0.5/1) and is a somatic-shaped signal - a diploid or haploid genotype model cannot express it; use a somatic caller (GATK Mutect2 --mitochondria-mode) for real heteroplasmy work. Polyploid/pooled samples need --ploidy N set to the true copy number so dosage/allele-count is preserved rather than collapsed to het.
A raw caller VCF is not a finished callset. The standard downstream order:
bcftools norm -f reference.fa -m -any variants.vcf.gz -Oz -o norm.vcf.gz. Do this before ANY comparison, annotation, or merge. See variant-calling/variant-normalization.QUAL, FORMAT/DP, SP) suited to the depth and platform. See variant-calling/filtering-best-practices.If the point of choosing bcftools vs a reassembly caller is accuracy, compare them correctly - this is where naive analyses go wrong:
bcftools norm -f ref.fa -m -any). Two VCFs can encode the identical haplotype with different records (indel placement in repeats, MNP vs split SNVs); un-normalized records mismatch spuriously.bcftools isec. A line-diff / isec on raw records overcounts both false positives and false negatives from representation alone. Score against a GIAB truth set with hap.py + vcfeval inside the confident-region BED (Krusche 2019), reporting SNVs and indels separately.Goal: Speed up calling on large inputs.
Approach: Pipe uncompressed BCF between stages, thread both tools, and shard by chromosome.
# Threaded, uncompressed-BCF pipe
bcftools mpileup -Ou -f reference.fa --threads 4 input.bam | \
bcftools call -mv --threads 4 -Oz -o variants.vcf.gz
# Parallel by chromosome, then concatenate
for chr in chr1 chr2 chr3; do
bcftools mpileup -Ou -f reference.fa -r "$chr" input.bam | \
bcftools call -mv -Oz -o "${chr}.vcf.gz" &
done
wait
bcftools concat -Oz -o all.vcf.gz chr*.vcf.gz
bcftools index all.vcf.gz
MQ <40 signal ambiguous mapping. -q 20 helps; a reassembly/alt-aware caller helps more.-d to 3-4x expected mean coverage; the default 250 truncates deep targeted panels and can bias likelihoods.| Symptom | Cause | Fix |
|---|---|---|
no FASTA reference |
-f omitted |
Add -f reference.fa |
[E::faidx] ... different number of sequences / reference mismatch |
mpileup reference != alignment reference | Use the exact FASTA the BAM was aligned to; compare @SQ in samtools view -H against grep '^>' ref.fa |
| No variants called | Coverage too low, -q/-Q too strict, empty/wrong BAM |
Check samtools depth; relax -q/-Q; confirm reference build |
| False heterozygous calls everywhere on chrX/chrY (male) | Non-PAR sex chromosome called as diploid | Set --ploidy 1 for non-PAR, or use a --ploidy-file / --ploidy GRCh38 |
| Excess indel false positives in repeats | Position-based limitation, not a bug | Normalize + hard-filter; validate or recall indels with a reassembly caller |
| Downstream tools disagree on the same variant | Records not normalized | bcftools norm -f ref.fa -m -any before comparing/merging/annotating |