How to Convert BAM to VCF: A Step-by-Step Guide for Variant Calling

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-11-11

By the end of this walkthrough you will have a normalized, filtered VCF file that lists the SNVs and small indels present in a BAM alignment, along with the commands and quality checks needed to trust that file. There is no single "convert" button. BAM to VCF is a variant calling step, and the tools, reference genome, and filtering choices all matter.

Before starting, have four things on hand: the aligned BAM, the matching reference FASTA with a .fai index and a sequence dictionary, a variant caller (bcftools or GATK), and roughly 8 GB of RAM and 4 CPU threads for a human whole-genome BAM (less for exomes or targeted panels). Reproducible pipelines such as ILIAD exist specifically to stitch together "alignment, variant calling, and/or build conversion" without needing to run each tool by hand [1], and you can model your manual workflow on the same logic. The commands below assume Linux with samtools, bcftools, and GATK 4 on the PATH.

Prerequisites and Environment

Every step in this guide depends on clean inputs. Confirm each item before touching a variant caller.

• samtools 1.16 or newer. Handles BAM sorting, indexing, and inspection.
• bcftools 1.16 or newer. Same version line as samtools, typically shipped together.
• GATK 4.2 or newer if you prefer HaplotypeCaller.
• A reference FASTA that matches the exact reference build the BAM was aligned against. Mixing GRCh37 with GRCh38 is the single most common reason a REF allele check fails later.
• Reference indices. You need ref.fa.fai (from samtools faidx) and ref.dict (from gatk CreateSequenceDictionary or samtools dict). GATK requires the dictionary. bcftools requires the .fai.
• Disk space. Bulky intermediates are real. A CRAM-to-BAM step or a MarkDuplicates step can each produce files roughly the size of the input BAM.

Verify the versions quickly:

samtools --version | head -1
bcftools --version | head -1
gatk --version | head -1

Expected output is one version string per tool, each starting with the tool name followed by a version number. A missing GATK version usually means Java is not configured properly, since GATK is a Java jar wrapper.

Understanding the File Formats First

Converting BAM to VCF only makes sense if you know what each format holds. A quick tour of the bioinformatics file formats involved prevents most of the mistakes covered later.

FASTQ vs FASTA

Both are plain text. The difference is quality data.

• FASTA stores a sequence per record, using a > header line followed by the sequence. There is no per-base quality. Reference genomes, transcripts, and protein databases are FASTA.
• FASTQ stores sequencing reads, with a four-line record per read: an @ identifier, the sequence, a + line, and an ASCII-encoded Phred quality string of the same length as the sequence.

A one-line illustration:

FASTA:  >chr1
        ACGTACGTAC...
FASTQ:  @read1
        ACGTACGT
        +
        IIIIFFFF

In the FASTQ-to-VCF data flow, FASTQ is the input to the aligner (BWA, Bowtie2, minimap2), BAM is the aligner output, and VCF is the caller output. If you start with FASTQ, you first align and sort before you can call variants.

SAM, BAM, and CRAM

SAM is the human-readable alignment format: tab-delimited lines with an eleven-field mandatory core plus optional tags. BAM is the BGZF-compressed binary form of SAM and is what most callers read. CRAM is a reference-compressed archival format that can shrink storage significantly, at the cost of needing the reference for full readout.

Large consortium projects have worked to make data access compatible across these formats. The Crypt4GH container, for example, adds byte-level random access to encrypted SAM, BAM, CRAM, VCF, and BCF files so that a tool can read sections without decrypting the whole file [2]. If you are handling human germline data under consent restrictions, that matters when you decide whether to keep BAM, CRAM, or both.

VCF and BCF

VCF is the tab-delimited plain-text variant format. Lines starting with ## are metadata, the single #CHROM line names the samples, and one row per variant follows. Each row has CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT, and one column per sample. BCF is the BGZF-compressed binary version of VCF, same information density at smaller file size and faster random access.

VCF is not the only downstream format. Pgltools, for instance, defines a separate paired-genomic-loci standard for Hi-C interactions because VCF does not express two-locus contact data naturally [3]. For SNV and indel calls, VCF is the universal interchange.

A minimal header of a real bcftools output looks like this:

##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=chr1,length=248956422>
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM  POS       ID  REF  ALT   QUAL  FILTER  INFO  FORMAT  SAMPLE
chr1    10177     .   A    AC    43.2  .       DP=9  GT:DP   0/1:9

Read that one row slowly. REF is A, ALT is AC (an insertion), the QUAL is 43.2, and the sample is heterozygous (0/1). Once you can parse a row like that, you can audit any VCF you generate.

The Whole Workflow at a Glance

The path from BAM to a usable VCF has seven stages:

1. Validate the BAM integrity.
2. Sort and index the BAM.
3. Mark duplicates (strongly recommended for PCR-based libraries).
4. Optionally recalibrate base qualities (GATK path only).
5. Run the variant caller to emit raw variant calls.
6. Normalize and split multiallelics.
7. Filter, validate, and compare to truth data.

Each stage has a concrete command and a concrete output check. Skipping stages is where most "why does my VCF look wrong" threads start.

Step-by-Step: Convert BAM to VCF

Step 1: Validate the BAM

Start by confirming the BAM is complete. A truncated BAM is a silent source of under-counted variants.

samtools quickcheck aligned.bam && echo "BAM OK"

quickcheck exits with status 0 when the file ends with a proper EOF marker. If it prints a warning such as [E::bam_hdr_read] EOF marker is absent. The input is probably truncated, do not proceed. Re-copy or re-download the BAM.

Also run:

samtools flagstat aligned.bam

Expected output is a block of counts, for example:

1000000 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
50000 + 0 duplicates
940000 + 0 mapped (94.00% : N/A)
960000 + 0 paired in sequencing
480000 + 0 read1
480000 + 0 read2
920000 + 0 properly paired (95.83% : N/A)

The mapped percentage and properly paired percentage tell you whether the alignment is reasonable. Any value in the low 80s or below for a modern human library is worth investigating before calling variants.

Step 2: Sort and Index

Variant callers require coordinate-sorted, indexed BAM. If your aligner already produced aligned.sorted.bam, skip to the index step. Otherwise:

samtools sort -@ 8 -m 2G -o aligned.sorted.bam aligned.bam
samtools index aligned.sorted.bam

The flags matter:

• -@ 8 uses 8 threads for compression and sorting.
• -m 2G gives each thread up to 2 GB of RAM for in-memory sorting. For very deep data, raising this reduces the number of temporary files.

samtools index creates aligned.sorted.bam.bai. Verify with:

samtools idxstats aligned.sorted.bam

You should see one line per contig in the header, with the length in column 2, then mapped and unmapped counts, plus a final * line for reads with no coordinates. If contigs are missing or the counts look empty for a contig you expected, the alignment step most likely used a different reference than the caller will.

Step 3: Mark Duplicates

PCR duplicates are identical fragments that inflate coverage and skew allele balance. Mark them before calling. GATK MarkDuplicates is the standard:

gatk MarkDuplicates \
  -I aligned.sorted.bam \
  -O aligned.md.bam \
  -M duplicates_metrics.txt
samtools index aligned.md.bam

sambamba markdup -t 8 aligned.sorted.bam aligned.md.bam is a faster alternative that produces the same flags. Picard MarkDuplicates uses positional arguments instead: picard MarkDuplicates I=... O=... M=....

Open duplicates_metrics.txt and look at the PERCENT_DUPLICATION line. Values under 5 percent suggest a clean library. Values over 20 percent suggest either very low input or an over-amplified library, both of which affect downstream sensitivity.

Step 4: Optional Base Quality Score Recalibration (GATK Only)

BQSR corrects systematic base-calling errors using known variant sites. It is optional and adds hours for a whole genome, but improves sensitivity for some sequencers.

You need known-sites VCFs, typically dbSNP, the Mills and 1000G gold-standard indels, and 1000G phased SNVs.

gatk BaseRecalibrator \
  -R reference.fa \
  -I aligned.md.bam \
  --known-sites dbsnp.vcf.gz \
  --known-sites mills_indels.vcf.gz \
  --known-sites 1000G_phase1.snps.high_confidence.vcf.gz \
  -O recal.table

gatk ApplyBQSR \
  -R reference.fa \
  -I aligned.md.bam \
  --bqsr-recal-file recal.table \
  -O aligned.recal.bam
samtools index aligned.recal.bam

If you do not have the known-sites VCFs, skip BQSR. Running BaseRecalibrator without them is not meaningful. bcftools does not use BQSR at all, and the mpileup model is robust to uncalibrated qualities so long as the caller thresholds are set conservatively.

Step 5a: Call Variants With bcftools

bcftools is fast, dependency-light, and outputs standards-compliant VCF. The two-stage pattern is mpileup | call.

bcftools mpileup \
  -f reference.fa \
  -a AD,DP,ADF,ADR,SP \
  -Ou \
  aligned.recal.bam | \
bcftools call \
  -mv \
  -Oz \
  -o sample.raw.vcf.gz

bcftools index -t sample.raw.vcf.gz

Flags explained:

• -f reference.fa supplies the reference. Required.
• -a AD,DP,ADF,ADR,SP adds allele depth, total depth, strand-specific depth, and phred-scaled strand bias to the output. Skip these only if you want a narrow file.
• -Ou writes uncompressed BCF to the pipe. Faster than text.
• -mv tells bcftools call to emit multiallelic calls and variants only. Drop -v to also record reference sites, which you almost never want for a raw call set.
• -Oz compresses the output as BGZF VCF.

Expected output header:

##fileformat=VCFv4.2
##bcftools_mpileupVersion=1.17
##bcftools_callCommand=call -mv -Oz -o sample.raw.vcf.gz

Then check the count:

bcftools view -H sample.raw.vcf.gz | wc -l

A human whole-genome at 30x will produce roughly 4 to 5 million raw variant records, more for high-coverage exomes. If you see zero records, read the error section below.

Step 5b: Call Variants With GATK HaplotypeCaller

HaplotypeCaller reassembles reads locally and is the reference implementation of the GATK Best Practices pipeline. For a single sample:

gatk HaplotypeCaller \
  -R reference.fa \
  -I aligned.recal.bam \
  -O sample.raw.vcf.gz

To produce a gVCF for later joint genotyping:

gatk HaplotypeCaller \
  -R reference.fa \
  -I aligned.recal.bam \
  -O sample.g.vcf.gz \
  -ERC GVCF

The -ERC GVCF flag switches to reference-confidence emitting, producing one record per block of non-variant positions in addition to variant records. The resulting file is typically three to five times larger than a direct VCF but allows cohort-level genotyping without storing BAMs.

For a cohort, gather gVCFs:

gatk GenomicsDBImport \
  --sample-name-map cohort.sample_map \
  --genomicsdb-workspace-path gendb \
  -L intervals.list

gatk GenotypeGVCFs \
  -R reference.fa \
  -V gendb://gendb \
  -O cohort.vcf.gz

The cohort.sample_map is a two-column tab file: sample name and gVCF path. intervals.list is a list of genomic regions, usually one line per chromosome or per scatter interval.

Step 6: Normalize and Split Multiallelics

Raw VCFs often represent the same indel multiple ways across samples. Normalization left-aligns and trims redundant bases, and splitting breaks ALT=G,T into two rows. Both are non-negotiable for downstream tools.

bcftools norm \
  -f reference.fa \
  -m -both \
  -Oz \
  -o sample.norm.vcf.gz \
  sample.raw.vcf.gz

bcftools index -t sample.norm.vcf.gz

The flags:

• -f reference.fa is required for left-alignment.
• -m -both splits multiallelics and joins biallelic representations. The opposite direction (-m +any) merges.
• -Oz compresses.

To check that every REF allele matches the reference, add -c e (error) instead of the default warn:

bcftools norm -c e -f reference.fa -m -both sample.raw.vcf.gz -Oz -o checked.vcf.gz

If -c e throws, you have a reference mismatch between the BAM and the FASTA. Stop and fix that before filtering, because every downstream filter is meaningless if REF is wrong.

Step 7: Filter and Validate

Raw calls contain false positives. Common filters:

bcftools filter \
  -e 'QUAL<20 || INFO/DP<10' \
  -s LowQual \
  -Oz \
  -o sample.filtered.vcf.gz \
  sample.norm.vcf.gz

bcftools index -t sample.filtered.vcf.gz

Two flags do the work:

• -e is an exclude expression. Variants that match it are flagged, not removed.
• -s LowQual names the soft-filter. Use -i (include) with an inverse expression if you prefer to keep only passing rows.

The || means OR. So this flags variants with QUAL under 20 or INFO DP under 10. Pass-only extraction is a separate command:

bcftools view -f PASS -Oz -o sample.pass.vcf.gz sample.filtered.vcf.gz

Variant counts at each stage tell you whether filtering is too aggressive:

bcftools stats sample.filtered.vcf.gz > sample.stats.txt
grep "^SN" sample.stats.txt

You will see lines like SN 0 number of samples: 1, SN 0 number of records: 4200000, and SN 1 number of SNPs: 3900000. Compare SN and non-SN counts. A SNP:indel ratio near 4:1 is typical for human. A ratio of 1:1 or worse suggests a caller or filter problem.

For a truth comparison in benchmarking, use bcftools isec:

bcftools isec -p compare -Oz truth.vcf.gz sample.pass.vcf.gz

Or use hap.py if you want precision, recall, and F1 against GIAB truth sets. Both approaches are standard.

Key Options at a Glance

Reading the Output: What to Expect

After a successful run, look at three things in the final VCF.

Header sanity. bcftools view -h final.vcf.gz | head -20 should show ##fileformat=VCFv4.x, the reference contigs, and the sample name. A missing ##contig line will break many downstream tools that try to fetch intervals.

Row count. bcftools view -H final.vcf.gz | wc -l gives a record count. Compare to the raw non-filtered count. If the drop is over 30 percent, your hard filters are probably too strict for the coverage depth.

Allele balance. bcftools query -f '[%AD]\n' final.vcf.gz | head gives per-sample allele depths for the first few variants. For a true heterozygous call at 30x, expect roughly 12 to 18 reads for each allele. Values like 2/30 in a heterozygous call signal either a duplicate problem or a variant caller threshold issue.

For a forensic check of one locus:

bcftools view -r chr1:10177 final.vcf.gz

That prints the single record at that position, or nothing if it was filtered out. Compare to samtools mpileup output for the same locus if you want to see the pileup that convinced the caller.

Common Errors and Fixes

Below are the messages you will actually hit, and the fix for each.

[E::bam_hdr_read] EOF marker is absent. The input is probably truncated The BAM is incomplete. Re-copy from source. Do not pipeline a truncated BAM through variant calling.

[E::hts_idx_load3] Could not load local index file 'aligned.bam.bai' The BAM was not indexed, or the .bai is stale because the BAM was re-sorted after indexing. Run samtools index aligned.bam again.

Could not retrieve index file for 'reference.fa' Missing the .fai. Run samtools faidx reference.fa.

A USER ERROR has occurred: Cannot read ... the sequence dictionary is missing or out of sync GATK needs ref.dict. Create it with gatk CreateSequenceDictionary -R reference.fa -O reference.dict.

[E::sam_parse1] CIGAR and query sequence are of different length The BAM is corrupted or was produced by an incompatible tool version. Re-align from FASTQ.

The REF allele at chrN:POS is not present in the reference The BAM was aligned to a different reference than the FASTA supplied to the caller, or the FASTA is not the same build. Verify with samtools view aligned.bam chrN:POS and compare to samtools faidx reference.fa chrN:POS-POS.

no records loaded or an empty VCF The most common cause is a chromosome naming mismatch: the BAM uses chr1 and the reference uses 1, or the reverse. bcftools view -H final.vcf.gz | head will show what contigs are actually in the output. If they are all empty, check bcftools mpileup with -r chr1 on a specific region to confirm a single-locus call works.

bcftools mpileup: unsupported type ... for INFO field Your bcftools version is older than the fields you requested in -a. Upgrade to a matched version of samtools and bcftools from the same release.

Common Mistakes and Limitations

Starting from an unsorted or stale BAM. A BAM that was re-sorted after the .bai was created will make a caller read the wrong blocks. Always re-index after any sort or reshuffle.

Using a reference different from the one used for alignment. Build mismatch is the most common silent error. Confirm the ##reference line in the BAM header with samtools view -H aligned.bam | grep @SQ, and match the contig names exactly.

Skipping duplicate marking on PCR libraries. Duplicates inflate apparent allele support and are a routine source of spurious somatic calls. Mark them on every PCR-based library. Duplicate marking is less critical for PCR-free libraries, but it costs little and it documents the library type.

Calling ONT or PacBio long reads with a short-read caller. Short-read callers expect short, paired reads and calibrated base qualities. Long-read callers such as sTELLeR exist specifically because short-read methods struggle with repetitive regions, and sTELLeR outputs VCF so it plugs into the same downstream formats [4]. Choose the caller to match the platform.

Treating a VCF as analysis-ready. A raw VCF is an intermediate. Filtering, normalization, and often false-positive review are part of the workflow, not an optional finishing touch.

Forgetting to normalize before intersecting with public datasets. Left-alignment differences produce false discordance against dbSNP or gnomAD. Normalize both sides.

Assuming a single VCF covers every variant class. Structural variants, copy number changes, mitochondrial heteroplasmy, and microsatellite expansions each use separate callers and separate downstream filters. The BAM to VCF path in this guide covers SNVs and small indels only.

Ignoring the workflow layer. Custom Bash wrappers are hard to reproduce. Modular systems exist precisely to keep the pieces organized. CAMP provides a standard module architecture for metagenomics pipelines so that tools can run independently or in sequence without breaking each other, and ILIAD does the same for the genotype-to-VCF path [5][1]. Even for a single BAM, writing your commands as a versioned script with pinned tool versions turns a fragile process into something another analyst can rerun.

Overlooking data handling when data are sensitive. Whole-genome BAMs and VCFs are not small. If the underlying data must be encrypted at rest, Crypt4GH enables direct read access to encrypted SAM, BAM, CRAM, VCF, and BCF without decrypting to a scratch disk, which addresses the leakage risk of decrypted copies [2]. That is not a variant-calling choice, but it is a file-format choice that many pipelines ignore until an audit.

Thresholds copied between projects. A QUAL threshold tuned for 30x short-read WGS will over-filter a 200x amplicon panel and under-filter a 10x exome. Re-tune with bcftools stats per dataset. There is no universal number.

Expecting the caller to catch every variant. Sensitivity is bounded by coverage, library complexity, and map quality. Poorly mapped regions, tandem repeats, and segmental duplications are undercalled by every short-read caller. Long-read orthogonal data is the usual remedy, not a different filter.

Frequently Asked Questions

Can you convert a BAM file to VCF without variant calling?

No. There is no metadata-only conversion. A BAM contains aligned reads and a VCF contains genotype calls, so the variant caller is the conversion step. samtools mpileup alone produces a pileup, not a VCF. You always need a statistical model that converts read evidence into variant calls and genotypes.

What is the difference between a BAM file and a VCF file?

A BAM is a compressed binary list of aligned reads, and a VCF is a text list of variant sites and genotypes. BAM is the input level, VCF is the summary level. A 30x whole-genome BAM is often 80-100 GB, and the matching VCF is 100-200 MB.

Should I use GATK or bcftools for BAM to VCF conversion?

Use GATK HaplotypeCaller if you are following GATK Best Practices or joining many samples through gVCFs. Use bcftools if you want fewer dependencies, faster runtime, and straightforward scripting. Both produce spec-compliant VCF. Choose one and keep it consistent within a project so you can compare call sets.

Can you call variants directly from FASTQ files?

Not directly. FASTQ must be aligned to a reference to produce a BAM first. The order is FASTQ to BAM through an aligner, then BAM to VCF through a caller. Some end-to-end pipelines hide the intermediate alignment, but the alignment still happens inside.

What is a gVCF and when should I use one?

A gVCF (genomic VCF) is a VCF that includes reference-confidence blocks in addition to variant calls. Use a gVCF when you plan to joint-genotype multiple samples with GATK. For a single-sample final call set, a plain VCF is smaller and simpler.

Why is my VCF empty or missing known variants?

The most common cause is a chromosome naming mismatch between the BAM and the reference, such as chr1 versus 1. The second most common is an unsorted or unindexed BAM. Check samtools view -H bam | grep @SQ and bcftools view -h vcf | grep ##contig to compare naming.

Do I need to filter a VCF after variant calling?

Yes, filtering is part of the workflow. Raw calls include low-quality sites, strand-biased calls, and sites with insufficient depth. A minimal filter on QUAL and INFO/DP removes most obvious noise. For clinical or publication work, use truth-set benchmarking with hap.py or bcftools isec to set thresholds empirically.

How does FASTQ vs FASTA affect the BAM to VCF workflow?

FASTA is the reference input and FASTQ is the read input. FASTA supplies the genome the caller compares against. FASTQ becomes BAM through alignment and only then becomes VCF. Mixing the two, or supplying the wrong FASTA build, is what produces REF allele errors at the calling stage.

Sources

1. ILIAD: a suite of automated Snakemake workflows for processing genomic data for downstream applications.
2. Crypt4GH: a file format standard enabling native access to encrypted data.
3. Pgltools: a genomic arithmetic tool suite for manipulation of Hi-C peak and other chromatin interaction data.
4. Detecting transposable elements in long-read genomes using sTELLeR.
5. CAMP: a modular metagenomics analysis system for integrated multistep data exploration.