10 Small Variant Calling: Concepts

Concepts

10.1 The problem

A pile of reads overlaps a position. Most of them say A, a few say G. Variant calling is deciding what that means: whether this sample differs from the reference here, what its genotype is, and how confident you are.

Framed that way, the difficulty is visible immediately. Every part of the evidence is unreliable. Bases are miscalled. Reads are misplaced — Section 6.5 covers why a confidently aligned read can still be in the wrong locus. The reference is one arbitrary haplotype among many, and a position where your sample differs from it is not thereby interesting. The caller has to weigh all of that and produce a discrete answer.

A variant call is not an observation. It is an inference about a genotype, conditional on a model of how the data were produced. Two callers disagreeing on a site are usually not disagreeing about the reads; they are disagreeing about the model, or about where to put the threshold. Treating a VCF as a record of what was seen, rather than as a set of conclusions someone drew, is the root of most misuse of one.

10.2 What a caller actually computes

Underneath the discrete call is something continuous, and it is the part worth understanding because it is what filtering later operates on.

For a diploid site the caller evaluates the possible genotypes — homozygous reference, heterozygous, homozygous alternate — and computes how likely the observed reads would be under each. Those genotype likelihoods are the real output. The genotype in the file is just whichever one won, and the quality attached to it is the margin by which it won.

Two numbers are routinely confused, and they answer different questions in the same way that alignment score and mapping quality do:

QUAL asks whether there is a variant here at all.
GQ asks, given that there is one, whether this is the right genotype.

A site can have overwhelming evidence of some variant and genuine ambiguity about which genotype it is — high QUAL, low GQ. The reverse happens too. Filtering on one when you meant the other is a common and silent error.

Why modern callers are haplotype-based

Early callers evaluated each position independently. That works until two variants are close enough to sit on the same read, where the correct question is which combinations of nearby variants explain the reads — because a sequence of changes that looks implausible position by position is often a single well-supported haplotype. Reassembling a window and comparing candidate haplotypes is why indel calling improved substantially, and it is what almost every current caller does.

10.3 Calling is the easy part; filtering is the hard part

A caller run with permissive settings will find nearly every real variant. It will also emit a very large number of false ones. The interesting work is separating them, and that step — not the calling — is where callsets differ most and where methodological choices hide.

The reason it is hard is that the errors are not random. They concentrate in homopolymers, in tandem repeats, in segmental duplications, near indels, and at the edges of capture targets. Anywhere the alignment was uncertain, the variants will be too. So a filter that works on the genome’s easy regions will be miscalibrated on the parts you were most likely to care about.

A filter is part of your method

Every threshold you set defines which regions of the genome your study can see. Filtering hard enough to make a callset look clean will remove real variation from repetitive and duplicated sequence — the same regions MAPQ filtering removes, for the same underlying reason. That is often the right trade. It is never a neutral one, and it belongs in the methods section rather than in a parameter file nobody reads.

10.4 Joint calling and the N+1 problem

Calling a cohort together is not the same as calling each sample separately and merging the results, and the difference is worth being precise about.

The advantage is evidence sharing. A site where one sample has weak evidence is interpreted differently if forty other samples are confidently variant there — that is a real signal, not a borderline artifact. Joint calling also distinguishes “this sample is homozygous reference here” from “this sample had no data here”, a distinction that vanishes when you merge single-sample files and which matters enormously for any downstream analysis that counts alleles.

The cost is combinatorial. If adding one sample requires recalling the entire cohort, a growing study becomes quadratic in effort. The standard resolution is to separate per-sample work from cohort work: each sample is processed once into an intermediate that retains the likelihoods at every position, and the joint step combines those intermediates. Adding a sample then costs one sample’s work plus one recombination, not a full recall.

10.5 The same variant, written several ways

This one is a representation problem, and it is responsible for a surprising share of apparent disagreement between callsets.

Consider a deletion of one base in a run of six identical bases. Which base was deleted? Every answer produces the same sequence, and every answer is a different VCF record. The same ambiguity arises with any indel adjacent to a repeat, and with variants that can be written either as one complex event or as several simple ones.

The field’s answer is normalization (Tan et al. 2015): a canonical form, in which a variant is left-aligned as far as it can shift, and reduced to the fewest bases that still express the change. Two normalized records for the same event are identical; two un-normalized ones may not be.

This is why comparing two VCFs is not a set operation on positions. Done naively, a caller that spelled an indel differently scores as one false positive plus one false negative, which is doubly wrong. Comparison has to happen at the level of the haplotype implied by the variants — do these two callsets assert the same sequence — rather than at the level of the records. Any benchmark that does not do this understates concordance, and understates it most in repetitive regions, where the spellings differ most.

10.6 What a benchmark actually proves

Callers are evaluated against reference materials — genomes characterized exhaustively enough to serve as an answer key (Zook et al. 2019). This is real progress, and it is routinely over-read.

The essential caveat is that a truth set comes with confident regions, and they are not the whole genome. They are the part where the truth is itself established with confidence, which by construction excludes much of what is hard. A caller reporting high accuracy has reported it on that subset. It has demonstrated very little about the regions that were excluded — and those are disproportionately the ones where callers disagree.

Two more things worth holding onto. Stratified results — by region type, by variant type — are far more informative than a single F-score, because performance is not uniform and the aggregate is dominated by the easy majority. And there is a circularity risk with trained callers: a model trained on material derived from a truth set, evaluated against that same truth set, will look better than it is.

The question a benchmark answers

Not “is this caller accurate” but “how did this caller do on this sample, in these regions, against a truth set assembled with these methods”. That is still a useful answer. It is just narrower than the number suggests, and the gap between the two is where most over-claiming lives.

10.7 What to take forward

A variant call is an inference about a genotype under a model, not a record of what was observed.
Genotype likelihoods are the real output; the call is the argmax and the quality is the margin. QUAL and GQ answer different questions.
Errors concentrate where alignment was already uncertain, so filtering is where callsets diverge and where methodological choices should be declared.
Joint calling shares evidence across samples and distinguishes homozygous reference from no data; doing it without recalling the cohort each time is the reason for per-sample intermediates.
One biological event has many valid spellings, so callsets must be normalized and compared by haplotype rather than by position.
A benchmark is a statement about confident regions, not about the genome.

Section 11.1 names what to run.