Bio Variant Calling

Version Compatibility

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

Variant Calling from a BAM

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

The Governing Principle

There are two families of short-variant caller, and the choice between them is the single most consequential decision here.

• Position-based genotype-likelihood callers (bcftools mpileup|call, the old samtools/UnifiedGenotyper lineage) trust the aligner's per-read placement. At each reference position they tally the pileup, compute P(reads | genotype) per site, and call under a Bayesian model. Fast, transparent, no training data. But the mapper places each read greedily and independently, so an indel near a read end or inside a repeat is placed inconsistently across reads, and the per-position model cannot repair that. bcftools mitigates it with BAQ (per-Base Alignment Quality: base qualities near a likely misalignment are downweighted so a shaky column does not produce a confident false SNP), but BAQ suppresses false positives rather than reconstructing the true indel.
• Local-reassembly / haplotype callers (GATK HaplotypeCaller, DeepVariant, DRAGEN) discard the local alignment in an active region and re-derive it: assemble candidate haplotypes, realign every read to them (PairHMM or a learned model), then genotype. The indel is represented once, on the assembled haplotype, instead of (mis)placed per read. This is precisely why they beat pileup callers on indels, clustered variants, and difficult regions.

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.

Engine Selection (the decision that comes before any command)

Guidance, not dogma; on a production human pipeline, validate against current GIAB/GA4GH benchmarks (hap.py + vcfeval) before committing.

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.

bcftools mpileup + call

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.

Basic calling

bcftools mpileup -f reference.fa input.bam | bcftools call -mv -Oz -o variants.vcf.gz
bcftools index variants.vcf.gz

Recommended single-sample pipeline

# -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

Region-restricted and multi-sample calling

# 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

The mpileup / call flags that change results

The multiallelic caller (-m) handles sites with several ALT alleles natively and is statistically superior; the consensus caller (-c) exists only for backward reproducibility.

Ploidy: sample, organelle, and sex-chromosome calling

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.