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.
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.
bcftools mpileup -f reference.fasta sorted.bam | bcftools call -mv -Oz -o calls.vcf.gzExpected output
[mpileup] 1 samples in 1 input files
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid2Index 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.
bcftools index calls.vcf.gz
bcftools view -H calls.vcf.gz | wc -lExpected output
48213Troubleshooting
[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`.