All protocols

Call variants with bcftools mpileup and call

Pile up read alignments and call SNPs and indels from a sorted BAM with bcftools, producing a compressed, indexed VCF.

SDSomenath DuttaUpdated July 21, 2026

Prerequisites

  • bcftools installed (`bcftools --version`)
  • The SAME reference FASTA used to align the reads
  • A coordinate-sorted, indexed BAM (sorted.bam plus sorted.bam.bai)

1Pile up and call in one pipe

Stream the pileup straight into the caller so no intermediate BCF touches disk. `-m` selects the multiallelic caller, `-v` keeps only variant sites, and `-Oz` writes a bgzipped VCF to calls.vcf.gz.

bash
bcftools mpileup -f reference.fasta sorted.bam | bcftools call -mv -Oz -o calls.vcf.gz

Expected output

[mpileup] 1 samples in 1 input files
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid

2Index the VCF and count the calls

Build a CSI index so downstream tools can random-access the file, then count records to sanity-check the run. `bcftools index` writes calls.vcf.gz.csi; the `view -H` pipe drops the header so `wc -l` counts only variant lines.

bash
bcftools index calls.vcf.gz
bcftools view -H calls.vcf.gz | wc -l

Expected output

48213

Troubleshooting

[E::faidx_fetch_seq] The sequence "chr1" was not found, or zero variants are called

The reference contig names must match the BAM header exactly. Use the same FASTA the reads were aligned to, not a differently-named build.

Wrong genotypes on a non-diploid organism

bcftools assumes diploid by default. Set the ploidy on the caller, e.g. `bcftools call -mv --ploidy 1` for a haploid genome, or pass a ploidy file with `--ploidy-file`.