Variant Calling Explained: Finding a Mutation in Sequencing Data

A variant call is a probabilistic statement about a genotype, produced from noisy reads whose alignment is itself uncertain. Reading the output correctly means knowing which of those uncertainties each field in the file is reporting.

From reads to a pileup

Sequencing produces short reads — commonly 150 bases — from random positions across the genome. Aligning them to a reference places each read at a coordinate, and the depth at a position is how many reads cover it. Thirty-fold mean coverage means the average base is covered by about thirty reads, though coverage is Poisson-ish and locally very uneven.

Stack the reads at one coordinate and you have a pileup. Suppose at chromosome 1, position 1,000,000, the reference base is C and thirty reads cover it: eighteen say C and twelve say T. The naive reading is a heterozygous C-to-T substitution with an allele fraction of 12 divided by 30, or 0.40.

Each of those base observations arrives with a quality score. Phred quality Q encodes an error probability as Q equals minus ten times the base-ten logarithm of the probability the base is wrong, so Q20 is one error in a hundred, Q30 one in a thousand, Q40 one in ten thousand. Each read also carries a mapping quality, on the same Phred scale, expressing the aligner’s confidence that the read came from this locus at all — which is a different and often larger source of error than the base call.

Why a count threshold fails

“Call a heterozygote if the alternate fraction is between 0.25 and 0.75” is the rule everyone writes first, and it is wrong in both directions for reasons that are easy to state.

• It ignores quality. Twelve alternate reads at Q40 and twelve at Q13 are wildly different evidence. At Q13 the per-base error rate is about 5 percent, and twelve errors out of thirty at that rate is far more plausible than twelve out of thirty at Q40.
• It ignores depth. Two reads out of four is a fraction of 0.50 and almost no evidence. Fifteen out of thirty is the same fraction and strong evidence. A rule on the ratio alone cannot distinguish them, which is precisely what a likelihood does.
• It assumes a diploid balance that is often absent.Somatic variants in a tumour sample sit at whatever fraction the clone occupies, which can be 0.05. Mosaic variants likewise. A heterozygote rule discards them by construction, which is why somatic callers are separate programs with different models rather than germline callers with a lower threshold.
• It ignores strand and position bias. Twelve alternate observations all on reads in the same orientation, or all within five bases of a read end, is a signature of an artefact rather than a variant. Counting cannot see it.

Genotype likelihoods, worked

The probabilistic formulation asks, for each possible genotype, how likely the observed bases are if that genotype were true. For a biallelic site with reference C and alternate T the three candidates are C/C, C/T and T/T, and the likelihood of one base observation given a genotype follows from the error probability e derived from its Phred score:

P(observe T | genotype C/C) = e / 3       // it is an error
P(observe T | genotype C/T) = 0.5 - e/3  // ~1/2, either allele
P(observe T | genotype T/T) = 1 - e      // it is correct

multiply across all 30 independent observations,
take logs to keep it numerically sane

With eighteen C and twelve T all at Q30, so e equal to 0.001, the C/T likelihood dominates by an enormous margin: producing twelve T observations under C/C requires twelve independent errors, each of probability about 0.00033, and producing eighteen C observations under T/T requires eighteen. The three likelihoods are then combined with a prior — human sites are mostly homozygous reference, and the prior differs between transitions and transversions and between known and novel sites — and normalised to posteriors.

Production callers add a step before this one, and it is the step that matters most for indels. Rather than assuming the input alignment is correct, they identify active regions showing evidence of variation, discard the existing alignment there, assemble the reads locally into candidate haplotypes with a de Bruijn graph, and realign every read against each candidate haplotype with a pair hidden Markov model. The reason is that an aligner placing a read against the reference has a strong bias toward calling a real insertion as a run of mismatches, because that is cheaper under its gap penalties. Local reassembly removes the reference from the middle of the decision.

What lands in the VCF

The Variant Call Format is a tab-delimited text file with eight fixed columns followed by per-sample columns. The specification is maintained by the GA4GH file-formats group and published as the VCF specification. One line for the site above:

#CHROM POS      ID  REF ALT QUAL FILTER INFO              FORMAT      SAMPLE1
chr1   1000000  .   C   T   612  PASS   DP=30;AF=0.40;MQ=60  GT:AD:DP:GQ:PL  0/1:18,12:30:99:612,0,540

• QUAL is the Phred-scaled probability that there is no variant at this site. 612 is an assertion about the site, not about the genotype.
• GT is the genotype: 0/1 heterozygous, 1/1 homozygous alternate, 0/0 homozygous reference, ./. no call. A vertical bar instead of a slash means the alleles are phased.
• AD is allelic depth as a comma-separated list, one entry per allele including the reference — here 18 reference reads and 12 alternate. DP is total depth, which is not always the sum of AD, because reads failing filters are counted in one and not the other.
• GQ is the Phred-scaled confidence in the called genotype specifically, capped at 99 by most callers. It is the field to filter on when you care about the genotype rather than the existence of the site, and it is routinely confused with QUAL.
• PL gives Phred-scaled likelihoods for each genotype in a defined order, normalised so the most likely is 0. Here 612,0,540 says C/T is best, C/C is 612 worse, T/T is 540 worse.

Filtering after calling is a separate stage with its own annotations: quality normalised by depth, strand-bias statistics, mapping quality, and the position of the alternate allele within reads. Setting hard thresholds on those is one approach; training a model on the annotations of known true sites and applying it to novel ones is the other.

Where calling is genuinely hard

Accuracy is not uniform across the genome, and the variation is structural rather than random. Short reads cannot be placed uniquely in segmental duplications and repeat arrays, so mapping quality collapses and calls there are unreliable regardless of depth. Indels in homopolymer runs are hard because the polymerase itself slips, so the error mode and the variant look identical. Regions of extreme GC content amplify poorly and end up under-covered. And any variant larger than a read is invisible to this approach entirely; structural variant detection uses read-pair orientation, split reads and coverage depth, and is a different problem.

The community handles this with benchmark sets: a reference sample with a high-confidence variant list and, critically, a set of regions where the truth is considered known. Reporting sensitivity and precision outside those regions is not meaningful, and comparing two callers only inside them systematically flatters both.

Where the computation stops

A variant call is a statement that a difference from the reference is present in the sample. It is not a statement that the difference causes anything. Deciding whether a variant is clinically significant is a separate, structured process carried out against published interpretation criteria by qualified people, drawing on population frequency, segregation, functional evidence and case data — and in a diagnostic setting it is performed in an accredited laboratory under a regulatory framework that governs the whole pipeline, including the software. Nothing on this page is a diagnostic method, and a caller output should not be presented to anyone as a result. If contamination is what is producing your unexpected calls, that is a separate diagnosis with its own signature; see detecting contamination in sequencing data.