Variant Calling Explained: From Aligned Reads to a VCF

A sorted BAM file is just a pile of aligned reads — by itself it makes no claims about where a sample differs from the reference genome. Variant calling is the process that turns that pile of reads into a confident list of SNPs and indels, written out as a VCF. This post walks through the conceptual pipeline connecting the two: how reads stack into a pileup, how a caller tells a true variant from a sequencing error, and what tools like GATK4 HaplotypeCaller, bcftools call, and DeepVariant are actually doing under the hood.

From Aligned Reads to a Pileup

Every variant caller starts from the same raw material: a coordinate-sorted, indexed BAM or CRAM produced by an aligner such as BWA or Bowtie2. The first conceptual step any caller performs, stated or not, is to build a pileup — for every genome position, collect every overlapping read and stack their bases, base qualities, and strand orientation into one column.

Picture a single genomic position covered by 40 reads. A pileup line for that position looks roughly like this:

chr7  140453136  A  40  ..,,.,,...,,..,,,g,,..,,,,.,,,,,.,,,,  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

Left to right: chromosome, position, reference base, depth (40), then one symbol per read — . or , for a match on the forward or reverse strand, a letter like g for a mismatch — followed by a base-quality string. Depth, base composition, and per-base quality are the three ingredients everything downstream is built on.

Real Variant, or Just a Sequencing Error?

Not every mismatch between a read and the reference is a variant. Illumina sequencing has a real base-call error rate — commonly around 1 in 1,000 bases at a Phred quality of 30 — and PCR, library prep, and alignment artifacts add more noise on top. A caller's central job is separating that background noise from a real, biological difference between the sample and the reference.

A handful of signals usually distinguish the two:

None of these signals is decisive alone. A single low-quality read supporting an alternate allele is almost certainly an error; forty reads split roughly 20/20 between reference and alternate, all decent quality, on both strands, is almost certainly a real heterozygous variant. Exactly where a caller draws the line between "certainly noise" and "certainly real" is what the statistical models inside every caller are built to do.

Allele Depth and Genotype Likelihoods

The single most useful number in this process is allele depth (AD): of the reads covering a position, how many support the reference allele versus each alternate allele. AD of 20,25 — 20 reads backing the reference, 25 backing an alternate — is the raw evidence a genotype call rests on. For a diploid organism, a ratio near 50:50 is exactly what you would expect from a true heterozygous site, since each chromosome copy contributes roughly half the reads. A ratio like 38:2 at the same total depth looks nothing like a real het call and points instead toward sequencing error, contamination, or a copy-number oddity.

Modern callers do not stop at eyeballing that ratio. They compute a genotype likelihood — the probability of observing this exact pileup under each possible genotype (homozygous reference, heterozygous, homozygous alternate) — using a statistical model that weighs base quality, mapping quality, and depth together. The genotype with the highest likelihood wins, and the gap to the next-best genotype becomes the genotype quality (GQ) reported in the VCF. Callers differ in how they get there: some compute genotype likelihoods more or less directly from the pileup, while GATK's HaplotypeCaller discards the pileup in a region and reconstructs candidate haplotypes from scratch before scoring genotypes against them — costlier, but noticeably more accurate around indels.

Germline vs Somatic Variant Calling

Before picking a caller, it helps to know what kind of variant you are looking for, because germline and somatic calling are statistically different problems.

• Germline variants are inherited or arise once in the germline, present in essentially every cell of the body. In a diploid sample, a real germline variant sits at an allele fraction near 0, 50, or 100 percent — no realistic in-between.
• Somatic variants arise in individual somatic cells, most commonly studied in cancer, where a tumor sample is a mixture of tumor and normal cells, often multiple subclones. A somatic variant can sit at any allele fraction from a few percent up to 100 percent, depending on tumor purity.

This matters because a caller tuned for germline variants assumes exactly two haplotypes per site, and will systematically miss or misclassify a mutation present in only 8 percent of reads — it looks too much like noise under a germline model. Somatic callers are built around that low-allele-fraction regime instead, typically by comparing a tumor sample directly against a matched normal to cancel out sequencing artifacts and shared germline variants alike.

Meet the Callers: HaplotypeCaller, bcftools call, and DeepVariant

With pileup-versus-error and germline-versus-somatic in place, the three tools bioinformaticians reach for most make a lot more sense.

GATK4 HaplotypeCaller

GATK4's HaplotypeCaller is the de facto standard for germline calling in research and clinical pipelines. Rather than call variants pileup-column by pileup-column, it flags "active regions" showing evidence of variation, locally reassembles the reads there into candidate haplotypes, and realigns each read against every candidate with a pair hidden Markov model to score how well it supports each one. That local reassembly is what makes it notably better than naive pileup callers at indels, since it is not thrown off by an indel's exact placement being ambiguous read by read — the tradeoff is speed. For cohorts, the usual pattern is to call each sample into a GVCF, then joint-genotype the whole cohort together, improving sensitivity for rare variants only weakly supported in any single sample.

bcftools call

bcftools's call command, paired with bcftools mpileup, is the more classical pileup-based approach: compute genotype likelihoods directly from the pileup using a Bayesian model, then call the most probable genotype. It skips local haplotype reassembly entirely, which makes it considerably faster and lighter than HaplotypeCaller, at some cost to accuracy in complex indel regions. For SNP-heavy work, quick surveys, non-human organisms without a mature GATK resource bundle, or resource-constrained environments, it remains a well-validated, reasonable choice.

DeepVariant

DeepVariant, from Google, takes a different approach entirely: it converts the pileup at each candidate site into an image-like tensor of read bases, qualities, and strand information, and hands that image to a convolutional neural network trained to classify the correct genotype — the same kind of architecture used for image recognition. It has consistently ranked near the top of precisionFDA and Genome in a Bottle truth-set benchmarks for germline SNP and indel accuracy. Because it is trained rather than hand-modeled, it also adapts to different sequencing platforms — Illumina, PacBio HiFi, Oxford Nanopore — simply by swapping in a platform-specific model, something statistical callers must be manually retuned for.

The conceptual difference above maps onto three quite different commands:

# GATK4: assemble haplotypes, then genotype (germline, one sample)
gatk HaplotypeCaller -R ref.fasta -I sample.bam -O sample.g.vcf.gz -ERC GVCF

# bcftools: pileup, then call genotype likelihoods, in one pipeline
bcftools mpileup -f ref.fasta sample.bam | bcftools call -mv -Oz -o sample.vcf.gz

# DeepVariant: pileup as an image, scored by a trained neural network
run_deepvariant --model_type=WGS --ref=ref.fasta --reads=sample.bam --output_vcf=sample.vcf.gz

From Pileup to VCF

Whichever caller produces them, genotype calls are written out the same way: one line per variant site, with the read evidence behind each sample's call packed into FORMAT fields like GT, AD, DP, and GQ. Everything covered above — pileups, allele depth, genotype likelihoods, germline versus somatic reasoning — exists to produce the numbers that land in those fields. Reading a VCF without knowing where those numbers came from is a bit like reading a lab report without knowing what the assay measured.

Practical Takeaways

• Every caller's real input is a pileup: stacked bases, qualities, and strand information at each genomic position, not the raw reads themselves.
• No single piece of evidence is decisive on its own. Real variants show consistent allele depth, decent base quality, and support from both strands; errors are usually isolated, low-quality, and one-sided.
• Germline and somatic calling are different statistical problems: germline calling assumes an allele fraction near 0, 50, or 100 percent, while somatic calling has to detect real variants at low, unpredictable allele fractions, typically against a matched normal sample.
• HaplotypeCaller trades speed for accuracy through local haplotype reassembly; bcftools call is faster and lighter through direct pileup-based genotype likelihoods; DeepVariant replaces both statistical models with a trained neural network and currently leads most germline benchmarks.
• The genotype fields you eventually read in a VCF — GT, AD, DP, GQ — are just the caller's final, compact summary of everything covered here.