How to Use bcftools mpileup for Variant Calling

The bcftools mpileup command works as the first half of a two-step variant calling pipeline: it reads aligned sequence data, stacks up all the reads covering each genomic position, and calculates genotype likelihoods. Those likelihoods then get piped into a second command, bcftools call, which makes the actual variant calls. The division of labor between these two tools is deliberate and worth understanding, because getting the best results means knowing what each step does and which flags matter at each stage.

The Two-Step Pipeline

Variant calling with bcftools is built around a pipe. You run bcftools mpileup to produce raw genotype likelihoods, and you pipe that output directly into bcftools call, which evaluates those likelihoods and decides which positions are genuine variants. The basic command looks like this:

bcftools mpileup -f reference.fa alignments.bam | bcftools call -mv -Oz -o output.vcf.gz

In this pipeline, bcftools mpileup reads your BAM file (or multiple BAM files) and, for every position in the genome, constructs a vertical “pileup” of all reads that overlap that position. It then calculates how consistent the observed bases are with each possible genotype, factoring in the mapping quality of each read, the base quality at each position, and the probability of local misalignment.¹ These per-position likelihood calculations are what bcftools call uses to decide whether a site is a variant.

The reason these live as two separate commands rather than one monolithic tool is flexibility. You can adjust the calling model, swap in different priors, or process multiple samples jointly, all by changing flags on bcftools call without touching the mpileup step. The pipe also keeps memory usage manageable, because the pileup data streams through rather than being written to a massive intermediate file.

What You Need Before You Start

Before running bcftools mpileup, you need two things: a reference genome in FASTA format and one or more alignment files in BAM or CRAM format. The reference genome must be the same one your reads were aligned to. It also needs to be indexed with samtools faidx, which creates a .fai index file. Your BAM files need to be sorted by coordinate and indexed as well (the .bai index files that samtools index produces).

The -f flag in the command above points to your reference FASTA. If you have multiple samples, you list all their BAM files after the reference, or provide them in a text file (one path per line) and pass that with -b. The mpileup step reads alignments from all files simultaneously, which is how it supports multi-sample calling.

Choosing a Calling Model

The bcftools call step offers two different statistical models, and the choice between them matters more than most people realize when they first set up their pipeline.

The original model, invoked with bcftools call -c, is a consensus caller that works well for simple biallelic sites (positions where there is one reference allele and one alternate allele). It was the default in the earliest versions of the samtools/bcftools variant calling pipeline.¹

The newer model, invoked with bcftools call -m, handles positions with multiple alternate alleles and also supports gVCF output.¹ For most current projects, -m is the better default. A benchmarking study comparing the multiallelic model to the consensus model found that the multiallelic approach produced more accurate calls in genetically diverse samples, such as those from African populations, where multiallelic sites are more common.² If your study organism or population has high genetic diversity, the multiallelic model is especially important.

The -v flag you see in most example commands tells bcftools call to output only variant sites, suppressing the thousands of reference-matching positions that would otherwise fill your VCF. If you want a record of every position, including non-variant ones (useful for distinguishing “no variant” from “no data”), drop the -v flag or use gVCF mode with -m.

Key Flags for bcftools mpileup

The mpileup step has a large set of options, but a handful matter most in practice.

-f (reference): Required. Points to your indexed FASTA reference genome.
-r (region): Restricts calling to a specific region, such as a single chromosome or a coordinate range like chr1:1000000-2000000. Useful for testing your pipeline on a small region before running the whole genome, and for parallelizing across chromosomes.
-q (minimum mapping quality): Filters out reads with mapping quality below the threshold you set. Reads that map ambiguously to multiple locations get low mapping quality scores, and including them introduces noise. A common starting point is -q 20.
-Q (minimum base quality): Filters out individual bases below a base-quality threshold. Sequencing errors are concentrated in low-quality bases, so raising this threshold reduces false positives at the cost of throwing away some data. A threshold of 20 is typical.
-d (max depth): Caps the number of reads considered at each position. The default is 250, which works for most whole-genome sequencing projects. If you have very high coverage data (targeted panels, amplicon sequencing), you may need to raise this.
-a (annotations): Tells mpileup which per-site annotations to compute and pass along in the output. Adding -a AD,DP gives you allelic depth and total depth fields, which are useful for downstream filtering.

These flags control the quality of the raw likelihoods that feed into bcftools call. Getting them right is more impactful than fiddling with post-hoc filters, because a bad pileup produces bad likelihoods that no amount of downstream filtering can fully rescue.

Base Alignment Quality and Why It Helps

One of the more subtle sources of false variant calls comes not from sequencing errors but from alignment artifacts. A read that is slightly misaligned at its edges can place the wrong base at a given position, and the base-quality score alone does not catch this because the sequencer reported the base correctly. The read just landed in the wrong spot.

Base Alignment Quality (BAQ) was developed specifically to address this problem. It measures the probability that each base in a read is wrongly aligned, and bcftools mpileup uses it to down-weight bases that are likely alignment artifacts.³ BAQ is enabled by default in bcftools mpileup, and in most cases you should leave it on. The main scenario where you might disable it (with -B) is when calling variants in reads that were aligned with a method already accounting for local realignment, or when working with data where the BAQ computation itself introduces problems, such as very short amplicon reads.

Filtering Your Call Set

Running the mpileup-to-call pipeline gives you a raw VCF, and that raw output will contain some false positives. Filtering is where you clean those up, and the approach with bcftools is more straightforward than with some competing tools.

The simplest and most effective first-pass filter is on the QUAL field, which represents the confidence that a variant is real. A benchmark study using simulated insect population data found that discarding SNVs with a QUAL score below 20 cut false positives by roughly 35 to 42 percent while dropping the recovery rate by less than one percentage point.⁴ That is a favorable trade-off. You can apply this filter with:

bcftools filter -i 'QUAL>=20' input.vcf.gz -Oz -o filtered.vcf.gz

Or you can use bcftools view with an include expression for the same effect. Beyond QUAL, you can filter on depth (removing sites with abnormally high or low coverage), strand bias (sites where the variant allele appears only on forward or reverse reads), and mapping quality. The bcftools filter and bcftools view commands support arbitrary expressions, so you can chain multiple conditions together. A reasonable starting set for whole-genome data might look like:

bcftools filter -i 'QUAL>=20 && DP>=5 && DP<=500 && MQ>=30' input.vcf.gz -Oz -o filtered.vcf.gz

The specific thresholds depend on your data: coverage depth, sequencing platform, and organism all affect what is reasonable. Treat these numbers as a starting point and adjust based on how your quality metrics respond.

Multi-Sample Calling

One of bcftools mpileup’s strengths is how naturally it handles multiple samples. You supply all your BAM files at once, and the pileup step considers evidence across all samples simultaneously when calculating likelihoods. The basic command extends to:

bcftools mpileup -f reference.fa sample1.bam sample2.bam sample3.bam | bcftools call -mv -Oz -o cohort.vcf.gz

Or, for many samples:

bcftools mpileup -f reference.fa -b bam_list.txt | bcftools call -mv -Oz -o cohort.vcf.gz

Joint calling across samples is important because it borrows strength: a variant that is weakly supported in one sample but clearly present in several others gets called more reliably than if each sample were processed alone. The bcftools call step uses allele frequencies estimated from the full set of input samples (or frequencies you supply externally) to evaluate genotypes under Hardy-Weinberg expectations.⁵

For very large cohorts (hundreds or thousands of samples), memory and runtime can become issues. A practical strategy is to split the genome into regions using -r, run mpileup | call on each region in parallel, and then concatenate the results with bcftools concat. This approach scales well on computing clusters and does not sacrifice calling accuracy as long as your region boundaries do not split individual variants.

How bcftools mpileup Compares to GATK HaplotypeCaller

GATK’s HaplotypeCaller is the other widely used variant caller, and people frequently ask which one to use. The two tools take fundamentally different approaches: bcftools mpileup reads the alignment as-is and builds pileup-based likelihoods, while HaplotypeCaller performs local reassembly of the reads into haplotypes before calling variants. The reassembly approach is powerful for complex regions but comes with higher computational cost and some quirks around filtering.

A benchmark study using simulated populations of a non-model insect found that bcftools mpileup had substantially lower false-positive rates than GATK HaplotypeCaller. The false-positive proportion for bcftools ranged from about 0.004 to 0.03 percent, compared to roughly 1.8 to 2.0 percent for GATK, a difference of 54 to 521 times. Recovery rates, meanwhile, were comparable between the two tools, both in the range of 89 to 90 percent.⁴

A particularly telling finding was that GATK’s variant quality scores did not cleanly separate true positives from false positives in this non-model organism context, making hard-filtering with GATK difficult. Bcftools mpileup’s QUAL scores, by contrast, provided a cleaner separation, meaning that a simple QUAL threshold was an effective filter.⁴ The authors concluded that bcftools mpileup may be the first choice for non-human studies.

That said, these results come from non-model organisms with less polished reference genomes, and GATK has been heavily optimized for human data with curated resource bundles. For human clinical genomics, GATK remains the standard in many pipelines. The honest advice is that both tools work well, and for non-human species or projects where you lack a well-annotated truth set for GATK’s variant quality score recalibration (VQSR), bcftools mpileup’s simpler filtering model is an advantage rather than a limitation.

Repetitive Regions and Common False-Positive Traps

Regardless of which variant caller you use, repetitive regions of the genome are the primary source of false-positive calls. In the benchmarking study referenced above, the vast majority of false positives from both bcftools mpileup and GATK HaplotypeCaller originated from repeats.⁶ Reads from repetitive regions map ambiguously, and even with mapping quality filters, some misplaced reads survive and generate spurious variant signals.

A practical mitigation is to exclude known repetitive regions from your call set using a BED file of repeat annotations and the -T ^repeats.bed flag in bcftools mpileup, or by filtering them out afterwards with bcftools view. If repeat annotations are not available for your organism, you can generate a rough set using tools like RepeatMasker or WindowMasker and apply those as a mask. The improvement in precision, especially for non-model organisms, can be dramatic.

Checking the Quality of Your Final Call Set

After filtering, you want some assurance that your call set is reasonable before feeding it into downstream analyses. One widely used sanity check for SNVs is the transition-to-transversion ratio (Ti/Tv). Transitions (changes between chemically similar bases, like A↔G or C↔T) happen more frequently in real biology than transversions (changes between dissimilar bases). For whole-genome sequencing of humans, the expected Ti/Tv ratio for true variants is around 2.0 to 2.1. A Ti/Tv ratio well below 2.0 suggests an excess of false positives, since random errors tend to produce transitions and transversions at more equal rates.

A study designing a variant quality control pipeline for whole-genome data found that applying successive filtering steps raised the Ti/Tv ratio from about 2.04 to roughly 2.16, with transversions being removed at a higher rate than transitions, consistent with preferential removal of false positives.⁷ Watching this ratio climb as you apply filters is a good sign that your filtering is cleaning up noise rather than randomly discarding real variants.

You can compute Ti/Tv easily with bcftools stats filtered.vcf.gz | grep "^SN" or by parsing the full stats output. Beyond Ti/Tv, bcftools stats gives you distributions of quality scores, depth, and other metrics that help you spot problems: an unusual spike in variants at very low depth, for instance, or an asymmetric distribution of variant quality scores that suggests systematic bias.

Working with Non-Human and Non-Model Organisms

A substantial share of bcftools mpileup users work outside human genomics, in ecology, agriculture, evolutionary biology, and entomology. The tool’s design actually suits these contexts well. It does not depend on a curated set of known variants for calibration the way GATK’s recommended best practices do, and its quality scores are interpretable enough that simple threshold filters work effectively even when you have no truth set to train a model against.

For non-model organisms, a few practical considerations apply. Reference genome quality matters more than you might expect: a fragmented assembly with many short contigs can produce artifacts at contig boundaries, where reads from the true genomic context get forced into an incomplete reference. If your reference is highly fragmented, consider excluding very short contigs or regions near contig ends from your call set. Ploidy is another factor: bcftools mpileup and bcftools call default to diploid, but many organisms are polyploid or have variable ploidy across their genomes. You can specify ploidy with --ploidy in bcftools call or provide a ploidy file with --ploidy-file for more complex cases like sex chromosomes or mixed-ploidy species.

Population-level studies in non-model organisms also benefit from joint calling across all individuals, as described above. The allele frequency estimates improve with more samples, and rare variants that would be missed in single-sample calling become visible when the evidence is pooled. For projects with dozens or hundreds of individuals of a non-model species, the bcftools mpileup | bcftools call pipeline with the -m flag is a reliable and computationally efficient default.

Putting the Pipeline Together

For a complete working example, suppose you have 20 insect samples aligned to a reference genome, and you want to call SNVs with reasonable quality. Your pipeline would look something like:

bcftools mpileup -f ref.fa -b bam_list.txt -q 20 -Q 20 -a AD,DP | bcftools call -m -v -Oz -o raw_variants.vcf.gz

bcftools filter -i 'QUAL>=20 && DP>=5' raw_variants.vcf.gz -Oz -o filtered_variants.vcf.gz

bcftools stats filtered_variants.vcf.gz > variant_stats.txt

The first line generates raw variant calls. The second applies basic quality and depth filters. The third gives you summary statistics including Ti/Tv to verify the call set looks reasonable. From there you can refine: exclude repetitive regions, adjust thresholds based on the stats output, or add strand-bias filters if your data shows asymmetry. The pipeline is intentionally modular. Each step reads standard VCF or BCF and writes standard VCF or BCF, so you can insert additional processing steps (normalization with bcftools norm, annotation with bcftools annotate) wherever they make sense for your project.

One practical note on output format: the -Oz flag produces compressed VCF (gzipped), which saves substantial disk space for whole-genome data. You will need to index the output with bcftools index before using it in downstream tools that require random access. The alternative -Ob flag produces BCF, a binary format that is faster for bcftools to read and write but less universally supported by other tools. For archiving and sharing, compressed VCF is the safe choice. For intermediate pipeline steps where only bcftools will touch the files, BCF can shave minutes off large runs.

1
PubMed Central. Twelve years of SAMtools and BCFtools
2
Briefings in Bioinformatics. Simulation of African and non-African low and high coverage whole genome sequence data to assess variant calling approaches
3
PubMed Central. Improving SNP discovery by base alignment quality
4
PubMed Central. The evaluation of Bcftools mpileup and GATK HaplotypeCaller for variant calling in non-human species
5
GigaScience. Twelve years of SAMtools and BCFtools
6
Scientific Reports. The evaluation of Bcftools mpileup and GATK HaplotypeCaller for variant calling in non-human species
7
Scientific Reports. Empirical design of a variant quality control pipeline for whole genome sequencing data using replicate discordance