Variant Calling Explained: How a BAM Becomes an Annotated VCF

Quick answer: Variant calling is the step that decides where a sample’s DNA differs from a reference and how much evidence supports each difference. The input is a coordinate-sorted, duplicate-marked BAM plus the same reference genome used for alignment; the output is a VCF, where each line is a claim about one position. Two tools cover most beginner work: GATK HaplotypeCaller for germline variants and bcftools mpileup for a lighter, scriptable alternative. The real craft is not running the caller but deciding which calls to believe, because an unfiltered VCF is a list of hypotheses, not results.

A companion to the bioinformatics file formats guide and the genome assembly guide: assembly builds the reference, alignment places the reads against it, and variant calling reads the differences.

The pipeline, in order

Variant calling is a short pipeline with three hard prerequisites, and skipping any of them produces calls you cannot trust:

• The reference must be identical to the one used for alignment. A mismatch of even one build (GRCh37 versus GRCh38, or a different contig set) shifts every coordinate silently. This is the single most common cause of a “correct” pipeline producing nonsense.
• The BAM must be coordinate-sorted, duplicate-marked, and indexed. Callers walk the genome in order and need random access.
• Read groups must be present. Base quality recalibration, joint genotyping, and downstream tools read the sample identity from @RG; missing read groups cause many tools to fail or to merge samples incorrectly.

From there the standard stages are:

• Mark duplicates (samtools markdup, gatk MarkDuplicates) so PCR copies are not counted as independent evidence.
• Recalibrate base qualities. Historically BQSR with a known-sites VCF. Current GATK guidance has moved away from it for most short-variant germline work, so treat this stage as version-dependent rather than mandatory.
• Call variants. Per-sample calling writes a gVCF (a VCF with reference-confidence blocks), which lets new samples be added to a cohort later without re-processing everything.
• Joint genotype across samples so the caller sees all samples at a site at once, which is what makes cohort genotypes comparable.
• Filter. Hard filters for single samples, VQSR or a trained model for large cohorts.

# Minimal germline path (GATK 4.x style)
gatk MarkDuplicates -I aln.bam -O dedup.bam -M metrics.txt
samtools index dedup.bam

gatk HaplotypeCaller -R ref.fa -I dedup.bam -O sample.g.vcf.gz -ERC GVCF
gatk GenomicsDBImport --sample-name-map cohort.map --intervals chr1:1-250000000 \
    --genomicsdb-workspace-path db
gatk GenotypeGVCFs -R ref.fa -V gendb://db -O cohort.vcf.gz

# Hard filtering for a single sample
gatk VariantFiltration -R ref.fa -V cohort.vcf.gz -O filtered.vcf.gz \
    --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0" --filter-name "LowQual"

What the VCF actually says

Every VCF line is one variant with three kinds of information, and knowing which is which is most of the format:

• Fixed columns describe the site: CHROM, POS, REF, ALT, QUAL, FILTER, INFO. REF and ALT are the alleles, so a line is literally “at this position the reference has this base and the sample has that one”.
• INFO carries site-level evidence: total depth (DP), allele balance (AF), strand bias metrics (FS, SOR), and mapping quality statistics (MQ). Filters are written here.
• FORMAT and sample columns carry per-sample genotype: GT is the genotype, AD the allele depths, DP the read depth at that sample, and GQ the confidence. Genotype 0/1 is heterozygous, 1/1 homozygous alternate, ./. missing.

The FILTER column is where honesty lives. PASS means the call survived your filters; anything else means it did not, and downstream tools ignore non-PASS records only if you tell them to. A pipeline that reports “12,000 variants” without stating the filter used has reported nothing.

Why unfiltered calls are mostly error

A raw caller emits a long tail of low-confidence records that come from sequencing error, mismapping in repetitive regions, paralogous reads collapsed together, and strand artefacts. Filters target exactly those signatures:

• Low QD (quality by depth) means the evidence is shallow relative to how confident the caller is, typical of a mismatch.
• High FS (Fisher strand) means the alternate allele appears mostly on one strand, which points at a sequencing or alignment artefact rather than a real variant.
• Low MQ means reads covering the site align ambiguously, so the site itself is suspect.
• Excess heterozygosity or high DP outliers flags collapsed duplications where two genomic copies are being read as one.

VQSR replaces thresholds with a trained classifier, but it needs a large, well-curated training cohort, and it is exactly the kind of step you should not use on three samples.

Germline, somatic, and structural variants are different problems

Beginners often assume one tool answers all variant questions. It does not:

Question	Typical approach	Key difference
Inherited (germline) SNVs and indels	HaplotypeCaller, DeepVariant, bcftools	Variant is present in most cells; genotype is the output
Tumour (somatic) variants	Mutect2 with a matched normal	Needs to subtract the patient’s own germline variants
Structural variants (large events)	Manta, Delly, Sniffles	Callers work on read pairs and split reads, not base-level pileups
Copy number changes	CNVkit, GATK CNV	Output is segments with ratios, not alleles

The consequence is practical: a somatic study that runs a germline caller will report the patient’s normal variants as tumour variants, which is a catastrophic but very common mistake.

Validation is the part people skip

No variant call is proven by the caller that produced it. Confidence comes from comparison against a truth set and an orthogonal method:

• Benchmarking against GIAB. The Genome in a Bottle reference samples have high-confidence variant truth sets, so you can run your whole pipeline on HG002 and report precision and recall instead of hoping.
• Tool agreement. Calls supported by two independent callers, or by both a variant caller and an orthogonal assay, carry meaningfully more weight.
• Manual inspection in a browser. Loading the BAM and VCF together in IGV shows whether the alternate allele actually appears in the reads at that position, which settles most arguments in seconds.

If you cannot validate, the honest phrasing is “calls were filtered with X and are reported as candidate variants” — not “we identified 12,000 mutations”.

Annotation is the next pipeline

A VCF with coordinates answers where; annotation answers what. Functional annotation (with tools such as VEP, snpEff, or ANNOVAR) adds gene, transcript, consequence, and population frequency information, and it is what turns a coordinate list into a sentence a biologist can read.

# Merge, filter and annotate
bcftools mpileup -f ref.fa dedup.bam | \
  bcftools call -mv -Oz -o calls.vcf.gz
bcftools filter -s LowQual -e 'QUAL<20 || DP<10' calls.vcf.gz | \
  bcftools view -f PASS -Oz -o filtered.vcf.gz
tabix -p vcf filtered.vcf.gz

Two annotation cautions matter more than the tool choice. First, the reference transcript set determines the consequence, so the same variant can be missense under one annotation release and synonymous under another. Second, population frequency filters are only as good as the population in the database: a variant that looks common in a reference panel may be genuinely rare elsewhere.

Common failure modes

• Wrong reference build. Every coordinate shifts. Check the contig naming (chr1 versus 1) and the exact assembly version before blaming the caller.
• Forgetting to mark duplicates. PCR duplicates inflate allele depth and manufacture false confidence.
• Confusing raw output with results. Reporting counts from an unfiltered VCF is the most common error in student projects.
• Copying hard-filter thresholds between organisms. Filter values tuned on human data do not transfer to bacteria or to non-model species.
• Ignoring the strand-bias and mapping-quality columns entirely. These two fields catch most artefacts, and skipping them is why some studies report indels in low-complexity regions.

Key takeaways

• Variant calling is a short pipeline with three hard prerequisites: matching reference, sorted and duplicate-marked BAM, and read groups.
• A VCF is a set of claims; the FILTER column and INFO metrics are where the confidence of each claim lives.
• Unfiltered calls are dominated by mismapping and error; filtering, or a trained model for large cohorts, is part of the method, not an optional cleanup.
• Germline, somatic, structural, and copy-number questions need different tools, and using the wrong class fails silently.
• Validate against a truth set such as GIAB and state the filters used, or describe your output as candidate variants rather than findings.

Challenge

Download a small publicly available BAM with a matched reference, run both bcftools and GATK on one chromosome, and compare the call sets. Then pick three variants found by only one caller, open them in IGV, and work out what the reads actually show. That comparison teaches more about variant calling than any pipeline diagram, because it puts the disagreement where it belongs: in the evidence.

Working on a variant analysis for a project or dissertation? Ampersand Academy offers one-to-one mentoring in bioinformatics and NGS analysis — part of Mahadhi, a development and digital marketing company.

FAQ