Calling Variants From Your Own WGS BAM With bcftools

By the end of this guide you will have a filtered, normalized and indexed VCF of single-nucleotide variants and short indels called from your own aligned whole-genome reads. You will also have a sense of how much to trust each call. You need three things: a coordinate-sorted BAM or CRAM file of your reads aligned to a reference, the exact reference FASTA those reads were aligned against, and bcftools 1.18 or newer along with samtools and htslib. A 30x human genome BAM is roughly 60–90 GB (CRAM is about a third of that), and you should budget 40 GB of scratch space and several hours of wall time. We assume a Linux box with at least 8 cores and 16 GB of RAM. Everything below is command-line work, and nothing in it is medical interpretation.

1. Understand what variant calling is doing

Variant calling is the step that converts a pile of aligned reads into a statement about your genome. At each position the caller asks a narrow question: given the bases and quality scores of the reads stacked here, and given the reference base, what genotype best explains the data? The output is a genotype (homozygous reference, heterozygous, homozygous alternate) plus a probability attached to that genotype. Everything downstream, from ancestry inference to pharmacogene annotation, inherits the errors made here.

bcftools splits this into two commands for a reason. bcftools mpileup walks the BAM and computes genotype likelihoods at every position with coverage, writing a per-site record with a PL field. bcftools call then applies a prior and a model across samples and decides which sites are variant. The separation matters because mpileup is the expensive I/O-bound part and call is cheap, so you can re-call with different priors without re-reading the BAM.

The model is pileup-based rather than assembly-based. bcftools evaluates columns of aligned bases and does not locally reassemble the reads. GATK HaplotypeCaller reassembles each active region into candidate haplotypes, which is why it handles indels and clustered variants better ¹. Benchmarking studies consistently find the two tools agree closely on high-confidence SNVs and diverge on indels and in repetitive regions ². Our view is that for SNVs on a clean 30x genome, bcftools is fast, transparent and good enough. For indels you care about, treat bcftools calls as a screen rather than an answer.

2. Prepare the reference and the BAM

The reference must be the same assembly and the same file used for alignment, including whether it has alt contigs and whether chromosomes are named chr1 or 1. A mismatch produces silent nonsense or an immediate contig error.

# index the reference (creates GRCh38.fa.fai)
samtools faidx GRCh38.fa

# confirm the BAM header references the same contigs
samtools view -H sample.bam | grep '^@SQ' | head -5

# index the BAM if needed
samtools index -@ 8 sample.bam

Check coverage before you spend hours calling. A genome that looks like 30x in the vendor report can be 22x after duplicate removal.

samtools coverage sample.bam | awk '$1 ~ /^chr[0-9]+$/ {print $1, $6, $7}'

Column 6 is percent of bases covered and column 7 is mean depth. If mean depth on autosomes is below about 20x, expect heterozygote dropout: a true het site with 8 reads has a meaningful chance of showing 8 reference reads and being called homozygous. Depth is the single variable that most strongly controls sensitivity ³.

If your data is CRAM, pass the reference explicitly with -T or set REF_PATH, otherwise decoding fails or silently stalls.

3. Run bcftools mpileup

This is the command we use for a single human sample at 30x:

bcftools mpileup \
  -f GRCh38.fa \
  -q 20 -Q 13 \
  -C 50 \
  -d 250 \
  -a AD,DP,SP,ADF,ADR \
  -Ou \
  sample.bam \
  > /dev/null   # placeholder; see piped version below

The flags earn their place. The -q 20 flag drops reads with mapping quality below 20, which removes most multi-mapping reads in segmental duplications. The -Q 13 flag drops individual bases with base quality below 13. The -C 50 flag applies the mapping-quality downgrade for reads carrying excessive mismatches, which suppresses false positives from misaligned reads. The -d 250 flag caps per-file depth so that a 3000x pileup in a collapsed repeat does not consume all your memory. The -a AD,DP,SP,ADF,ADR flag requests allelic depths, total depth, strand-bias p-value, and forward and reverse allelic depths, all of which you need for filtering later.

-Ou writes uncompressed BCF to stdout, which is the right output when you are piping into bcftools call. Compressing and then immediately decompressing wastes cycles.

Do not run this on the whole genome as one process. Split by region and run in parallel:

# build a region list of autosomes plus X, Y, M
cut -f1,2 GRCh38.fa.fai | grep -P '^chr([0-9]+|X|Y|M)\t' \
  | awk '{print $1":1-"$2}' > regions.txt

call_region() {
  r=$1
  bcftools mpileup -f GRCh38.fa -r "$r" \
      -q 20 -Q 13 -C 50 -d 250 -a AD,DP,SP,ADF,ADR -Ou sample.bam \
  | bcftools call -m -v -Oz -o "raw.${r//:/_}.vcf.gz"
  bcftools index -t "raw.${r//:/_}.vcf.gz"
}
export -f call_region
parallel -j 8 call_region :::: regions.txt

On 8 cores this finishes a 30x genome in roughly two to four hours, dominated by chromosome 1 and 2. The --threads flag in bcftools only parallelizes compression, not the pileup itself, so region splitting is how you get real speedup.

4. Run bcftools call and pick a model

bcftools call has two models. -m is the multiallelic and rare-variant caller and is the one to use. -c is the original consensus caller from samtools, kept for backward compatibility, and it handles multiallelic sites poorly. Use -m.

bcftools call -m -v -Oz -o raw.vcf.gz

-v outputs variant sites only. Leave it off if you want all sites, which you do need for some downstream purposes. A gVCF-style all-sites file lets you distinguish “reference at this position” from “no data at this position”, and that distinction matters when you are checking whether a specific pharmacogenomic or carrier-screening position was covered at all. For an all-sites run use bcftools mpileup without -v on call and add -g 5 to emit compressed reference blocks:

bcftools call -m -g 5 -Oz -o allsites.g.vcf.gz

One prior worth knowing about: -P sets the expected mutation rate, default 1.1e-3. Lowering it makes the caller more conservative and raises it makes the caller more permissive. We leave it at the default for a human genome against GRCh38. Changing it is a blunt instrument compared with filtering on the annotations you already collected.

If you produced per-region files, concatenate them:

ls raw.chr*.vcf.gz | sort -V > files.txt
bcftools concat -f files.txt -a -Oz -o raw.vcf.gz
bcftools index -t raw.vcf.gz

5. Normalize, then filter

Normalization comes before filtering, because unnormalized indels compare and annotate incorrectly. Left-align indels against the reference and split multiallelic records into biallelic ones:

bcftools norm -f GRCh38.fa -m -both -Oz -o norm.vcf.gz raw.vcf.gz
bcftools index -t norm.vcf.gz

-m -both splits both SNVs and indels at multiallelic sites. -f enables left-alignment. Watch the stderr summary: it reports how many records were realigned and how many were skipped because the reference allele did not match, which is your check that you used the right FASTA.

Now filter. bcftools has no machine-learned filter equivalent to GATK’s variant quality score recalibration, so you write hard filters. These are the thresholds we use as a starting point for a single 30x germline genome:

bcftools filter -s LowQual -e 'QUAL<30' -m + norm.vcf.gz -Ou \
| bcftools filter -s LowDepth -e 'FORMAT/DP<10' -m + -Ou \
| bcftools filter -s HighDepth -e 'FORMAT/DP>200' -m + -Ou \
| bcftools filter -s StrandBias -e 'FORMAT/SP>30' -m + -Ou \
| bcftools filter -s LowMQ -e 'MQ<40' -m + -Oz -o filtered.vcf.gz
bcftools index -t filtered.vcf.gz

-s labels rather than deletes, and -m + appends labels instead of replacing them, so a site can carry several reasons. Keeping the failures in the file is better than dropping them: when a downstream annotation surprises you, you want to see whether the underlying call was marginal.

Indels deserve a stricter threshold than SNVs. A separate pass is clearer than one long expression:

bcftools filter -i 'TYPE="indel"' filtered.vcf.gz -Ou \
| bcftools filter -s IndelLowQual -e 'QUAL<50 || IDV<5' -m + -Oz -o indels.vcf.gz

Count what survived:

bcftools view -f PASS filtered.vcf.gz | bcftools stats - | head -30

For a single human genome against GRCh38 expect roughly 4.0 to 5.0 million passing SNVs and 400,000 to 700,000 short indels, with the exact number depending on ancestry and on how aggressively you filtered. The ratio of transitions to transversions should land near 2.0 to 2.1 genome-wide. A ratio below about 1.8 suggests you are keeping a population of false positives.

6. Check the calls against a truth set

You cannot evaluate a pipeline on your own genome without a reference answer, so evaluate it on one where the answer is known. The Genome in a Bottle consortium publishes high-confidence variant calls and confident-region BED files for several reference samples, most notably HG002 ⁴. Run your exact pipeline on a public HG002 BAM subset, then compare:

bcftools isec -p isec_out -n=2 -c none \
  filtered.hg002.vcf.gz HG002_GRCh38_benchmark.vcf.gz

For a real precision and recall number use a tool that handles representation differences, such as hap.py or vcfeval. Raw isec will count a left-aligned and a right-aligned representation of the same indel as two different variants. Benchmarking work has repeatedly shown that apparent discordance between callers is often representation rather than genuine disagreement, and that the hard part of benchmarking now lies outside the easy regions ⁴.

Expect SNV precision and recall above 0.99 inside the confident regions and a clear drop for indels, particularly insertions longer than about 10 bp. Comparisons across pipelines find that concordance is high in easy regions and that disagreement concentrates in repeats, homopolymers, and segmental duplications ²⁵. If you need trustworthy indel calls in a specific gene, use an assembly-based caller there and compare.

Two limits are worth stating plainly. First, this pipeline calls germline SNVs and short indels only. Structural variants, copy-number changes, repeat expansions, and mosaic or somatic variation need different tools and different statistics ⁶⁷. Second, a variant call is a measurement. It is not an interpretation. Any call you intend to act on should be confirmed by a clinical-grade laboratory and discussed with a clinician or genetic counselor.

7. BCF versus VCF, and reading the file

VCF is the tab-delimited text format. BCF is its binary equivalent with the same information model. They convert losslessly in both directions. BCF parses far faster because fields are typed and fixed-width and strings are dictionary-encoded, so a bcftools query over a BCF can be several times faster than over the same data as bgzipped VCF. Use uncompressed BCF (-Ou) between piped commands, compressed BCF (-Ob) for files you will query repeatedly, and bgzipped VCF (-Oz) for files you will share or feed to tools that only read text.

bcftools view -Ob -o filtered.bcf filtered.vcf.gz   # VCF -> BCF
bcftools index filtered.bcf                          # .csi index
bcftools view -Oz -o out.vcf.gz filtered.bcf         # BCF -> VCF

To open and inspect a VCF, never use a plain text editor on a multi-gigabyte file. Use bcftools:

bcftools view -h filtered.vcf.gz | less          # header only
bcftools view -H filtered.vcf.gz | head          # records only
bcftools view filtered.vcf.gz chr7:117,480,000-117,670,000   # one region
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t[%GT\t%AD]\n' \
  -r chr1:1000000-1010000 filtered.vcf.gz

bcftools query with a format string is the fastest route into a data frame. Pipe it to a TSV and load with pandas or Arrow rather than parsing VCF in Python.

Common problems

The reference allele does not match. bcftools norm reports a count of records where the REF field disagrees with the FASTA. Any nonzero count means you used a different reference than the aligner did, or a different patch release. Stop and fix the reference. Do not filter around it.

Contig naming mismatch. chr1 versus 1, or a BAM with alt contigs against a primary-only FASTA. bcftools errors out at startup for a full mismatch but a partial one can pass and quietly drop regions. Compare samtools view -H against the .fai.

Too many heterozygous calls with skewed allelic depth. If AD frequently reads something like 40,4 at het calls, you are picking up alignment artifacts in paralogous regions. Raise -q to 30 and check that the strand-bias filter is doing work.

Too few heterozygous calls. Heterozygosity below about 0.0008 per base on autosomes usually means low depth or an over-strict depth filter. Lower FORMAT/DP to 8 and re-examine. This is a sensitivity problem, and calling parameters have a measurable effect on downstream relatedness and genotype-concordance estimates ⁸.

Indels look wrong in homopolymers. Pileup-based callers systematically struggle in homopolymer runs and short tandem repeats. Intersect your calls with a low-complexity BED file and treat anything inside it as provisional ⁹.

Wrong assembly for a downstream tool. If an annotation database expects GRCh37 and your calls are GRCh38, convert with bcftools +liftover rather than a naive coordinate map, because it handles reference-allele swaps and strand flips that a simple chain lift gets wrong ¹⁰.

Memory blowup on mpileup. Almost always a missing or too-high -d cap hitting a pileup in centromeric or mitochondrial sequence. Keep -d 250 for a 30x genome and raise it only for targeted deep sequencing.

Footnotes

1. Geraldine A. Van der Auwera, Mauricio O. Carneiro, Christopher Hartl, et al. From FastQ Data to High‐Confidence Variant Calls: The Genome Analysis Toolkit Best Practices Pipeline. Current Protocols in Bioinformatics, 2013. https://doi.org/10.1002/0471250953.bi1110s43 ↩

2. Chengyue Tan, Weifen Sun, Jiaxin Yu, et al. Comprehensive assessment of whole-genome sequencing-based variant calling pipelines in monozygotic twins. Briefings in Bioinformatics, 2025. https://doi.org/10.1093/bib/bbaf652 ↩ ↩²

3. Stepanka Zverinova, Victor Guryev. Variant calling: Considerations, practices, and developments. Human Mutation, 2021. https://doi.org/10.1002/humu.24311 ↩

4. Nathan D. Olson, Justin Wagner, Nathan Dwarshuis, et al. Variant calling and benchmarking in an era of complete human genome sequences. Nature Reviews Genetics, 2023. https://doi.org/10.1038/s41576-023-00590-0 ↩ ↩²

5. Kirti Bhadhadhara, Muthukumar Balamurugan, Neeraj Bharti, et al. “Performance Evaluation of Variant Calling Tools for Human and Microbial Genomes”. 2023 International Conference on Emerging Trends in Networks and Computer Communications (ETNCC), 2023. https://doi.org/10.1109/etncc59188.2023.10284935 ↩

6. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩

7. Chang Xu. A review of somatic single nucleotide variant calling algorithms for next-generation sequencing data. Computational and Structural Biotechnology Journal, 2018. https://doi.org/10.1016/j.csbj.2018.01.003 ↩

8. August E. Woerner, Sammed Mandape, Kapema Bupe Kapema, et al. Optimized variant calling for estimating kinship. Forensic Science International: Genetics, 2022. https://doi.org/10.1016/j.fsigen.2022.102785 ↩

9. Qianqian Song, Taobo Hu, Baosheng Liang, et al. cascAGS: Comparative Analysis of SNP Calling Methods for Human Genome Data in the Absence of Gold Standard. Interdisciplinary Sciences: Computational Life Sciences, 2024. https://doi.org/10.1007/s12539-024-00653-8 ↩

10. Giulio Genovese, Nicole B Rockweiler, Bryan R Gorman, et al. BCFtools/liftover: an accurate and comprehensive tool to convert genetic variants across genome assemblies. Bioinformatics, 2024. https://doi.org/10.1093/bioinformatics/btae038 ↩