Bio Variant Calling

bcftools mpileup/call uses a Bayesian statistical model on pileup data. It is fast and effective for germline SNPs, but its indel calling is less accurate than local-assembly callers (GATK HaplotypeCaller, DeepVariant) in homopolymer and low-complexity regions. For production human germline pipelines, prefer DeepVariant or GATK; use bcftools for rapid exploration, non-model organisms, or when computational resources are limited.

Basic Workflow

BAM file + Reference FASTA
         |
         v
   bcftools mpileup (generate pileup)
         |
         v
   bcftools call (call variants)
         |
         v
   VCF file

bcftools mpileup + call

Goal: Detect SNPs and indels from aligned reads using the bcftools pileup-and-call pipeline.

Approach: Generate per-position pileup likelihoods with mpileup, then call genotypes with the multiallelic caller.

"Call variants from my BAM file" → Generate genotype likelihoods from aligned reads and identify variant sites using a Bayesian caller.

Basic Variant Calling

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

Output Compressed VCF

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

Call Specific Region

bcftools mpileup -f reference.fa -r chr1:1000000-2000000 input.bam | \
    bcftools call -mv -o region.vcf

Call from Multiple BAMs

bcftools mpileup -f reference.fa sample1.bam sample2.bam sample3.bam | \
    bcftools call -mv -o variants.vcf

BAM List File

# bams.txt: one BAM path per line
bcftools mpileup -f reference.fa -b bams.txt | bcftools call -mv -o variants.vcf

mpileup Options

Goal: Control pileup generation with quality thresholds, annotations, and region restrictions.

Approach: Set minimum mapping/base quality, request specific FORMAT/INFO tags, and restrict to target regions.

Quality Filtering

bcftools mpileup -f reference.fa \
    -q 20 \           # Min mapping quality
    -Q 20 \           # Min base quality
    input.bam | bcftools call -mv -o variants.vcf

Annotate with Read Depth

bcftools mpileup -f reference.fa -a DP,AD input.bam | bcftools call -mv -o variants.vcf

Full Annotation Set

bcftools mpileup -f reference.fa \
    -a FORMAT/DP,FORMAT/AD,FORMAT/ADF,FORMAT/ADR,INFO/AD \
    input.bam | bcftools call -mv -o variants.vcf

Target Regions (BED)

bcftools mpileup -f reference.fa -R targets.bed input.bam | \
    bcftools call -mv -o variants.vcf

Max Depth

bcftools mpileup -f reference.fa -d 1000 input.bam | bcftools call -mv -o variants.vcf

call Options

Calling Models

The multiallelic caller (-m) handles sites with multiple ALT alleles natively and is statistically superior. The consensus caller (-c) is retained only for reproducibility of legacy analyses.

Output Variants Only

bcftools mpileup -f reference.fa input.bam | bcftools call -mv -o variants.vcf
# -v outputs variant sites only (not reference calls)

Output All Sites

bcftools mpileup -f reference.fa input.bam | bcftools call -m -o all_sites.vcf
# Without -v, outputs all sites including reference

Ploidy

Correct ploidy is critical -- calling a diploid organism as haploid halves heterozygous sensitivity; calling haploid as diploid produces false heterozygous calls.

# Haploid calling (bacteria, mitochondria, chrY in males)
bcftools mpileup -f reference.fa input.bam | bcftools call -m --ploidy 1 -o variants.vcf

# Ploidy file for mixed-ploidy organisms (e.g., sex chromosomes, polyploids)
# Format: CHROM FROM TO SEX PLOIDY
bcftools mpileup -f reference.fa input.bam | bcftools call -m --ploidy-file ploidy.txt -o variants.vcf

# Built-in ploidy presets
bcftools call -m --ploidy GRCh38 ...  # Human with sex chromosome handling

Prior Probability

# Adjust variant prior (default 1.1e-3 for human)
# Lower for highly inbred lines; higher for diverse/outbred populations
bcftools mpileup -f reference.fa input.bam | bcftools call -m -P 0.001 -o variants.vcf

Common Pipelines

Goal: Run production-ready variant calling workflows for single-sample and multi-sample analyses.

Approach: Chain mpileup and call with quality filters, annotations, and compressed output, optionally parallelized by chromosome.

Standard SNP/Indel Calling

bcftools mpileup -Ou -f reference.fa \
    -q 20 -Q 20 \
    -a FORMAT/DP,FORMAT/AD \
    input.bam | \
bcftools call -mv -Oz -o variants.vcf.gz

bcftools index variants.vcf.gz

Multi-sample Calling

bcftools mpileup -Ou -f reference.fa \
    -a FORMAT/DP,FORMAT/AD \
    sample1.bam sample2.bam sample3.bam | \
bcftools call -mv -Oz -o cohort.vcf.gz

bcftools index cohort.vcf.gz

Calling with Regions

bcftools mpileup -Ou -f reference.fa \
    -R targets.bed \
    -a FORMAT/DP,FORMAT/AD \
    input.bam | \
bcftools call -mv -Oz -o targets.vcf.gz

Parallel by Chromosome

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

# Concatenate results
bcftools concat -Oz -o all.vcf.gz chr*.vcf.gz
bcftools index all.vcf.gz

Annotation Tags

INFO Tags

FORMAT Tags

Request Specific Annotations

bcftools mpileup -f reference.fa \
    -a FORMAT/DP,FORMAT/AD,FORMAT/SP,INFO/AD \
    input.bam | bcftools call -mv -o variants.vcf

Performance Options

Goal: Speed up variant calling for large datasets.

Approach: Use multi-threading and uncompressed BCF piping to reduce I/O overhead.

Multi-threading

bcftools mpileup -f reference.fa --threads 4 input.bam | \
    bcftools call -mv --threads 4 -o variants.vcf

Uncompressed BCF for Speed

bcftools mpileup -Ou -f reference.fa input.bam | bcftools call -mv -Ou | \
    bcftools filter -Oz -o filtered.vcf.gz

Quick Reference

Difficult Region Considerations

bcftools pileup-based calling is particularly vulnerable in:

• Homopolymers: Runs of identical bases cause indel false positives. Validate indel calls in homopolymers >6bp with a second caller or manual review.
• Low-complexity regions: ~35 Mb of GRCh38 is flagged as LCR. Consider excluding with -T ^lcr.bed or flagging calls in these regions.
• Segmental duplications: Reads from paralogs pile up, creating false SNPs. Check MQ values -- sites with mean MQ <40 suggest ambiguous mapping.
• High-depth regions: Set -d (max depth) to 3-4x the expected mean coverage. The default (250) may be too low for targeted panels or too high for low-coverage WGS.

Common Errors

Related Skills

• vcf-basics - View and query resulting VCF
• filtering-best-practices - Filter variants by quality
• variant-normalization - Normalize indels
• variant-calling/gatk-variant-calling - Higher accuracy with local assembly
• variant-calling/deepvariant - Deep learning alternative
• alignment-files/pileup-generation - Alternative pileup generation