11 Small Variant Calling: Implementations
Everything on this page is generated from tools/variants.json. If something here is wrong or out of date, the fix belongs in that file, not in this prose.
For what a variant call is and why filtering is the hard part, see Section 10.1.
11.1 Start here
The choice is driven by the data, not by a ranking. Short-read human cohorts, long reads, and everything that is neither have different answers, and the tools below are the defaults for each rather than competitors for one crown.
What samtools is to alignments, this is to variants: reading, filtering, annotating, merging and subsetting VCF and BCF, released in lockstep with samtools and htslib. Its filtering expression language is the part worth learning first, because it replaces a great deal of fragile text processing. (Danecek et al. 2021, 2011)
Built for long reads specifically, combining a fast pileup pass with a full-alignment pass on the sites the first pass finds ambiguous. That two-stage design is the reason it stays tractable on ONT data, where the error profile makes every site look ambiguous to a single-stage caller. (Zheng et al. 2022)
Replaces the hand-built statistical model with a classifier trained on images of the pileup, which means its notion of a convincing variant is learned from truth sets rather than specified. The practical consequence is that it ships separate models per data type — short reads, PacBio, ONT — and using the wrong one is the main way to get poor results from it. (Poplin et al. 2018)
The reference implementation of joint calling, and the pipeline most published human work is specified against. HaplotypeCaller emits per-sample intermediates that are combined across a cohort, which is what lets a new sample be added without recalling everyone — the property that matters once a study has more than a handful of genomes. (McKenna et al. 2010)
Not an aligner, but every alignment step ends in it: sorting, indexing and converting SAM to BAM or CRAM. In practice the aligner and the sort are one piped command, so the unsorted intermediate never reaches disk. (Danecek et al. 2021)
11.2 The other callers
Reached for when the defaults do not fit: a non-diploid organism, a species with no trained model, a design the mainstream tools do not cover, or a pipeline that already specifies one of these. Two of the three below are no longer moving, which is stated in each record rather than implied.
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| FreeBayes | active | C++ |
Haplotype-based and genuinely simple to run: one command, no model files, no per-platform training. That makes it a reasonable default for non-human organisms and for ploidies other than two, where the trained callers have no matching model and the joint-calling machinery is overkill. (Garrison and Marth 2012) |
| Octopus | dormant | C++ |
Calls germline, somatic and de novo variants through one haplotype model rather than three separate tools, which is a genuinely different design and an attractive one for trio and tumor-normal work. Listed as viable rather than recommended because of its release history, not its method. (Cooke et al. 2021) |
| Strelka2 | deprecated | C++ |
Fast and accurate for its era, and still embedded in a great many production pipelines, particularly for tumor-normal calling. It is here because you will meet it, not because a new pipeline should adopt it: the repository is archived upstream and has had no release since 2018. (Kim et al. 2018) |
11.3 Evaluating a callset
Section 10.5 is the reason this section exists as something other than a footnote. Comparing callsets by position rather than by haplotype produces numbers that are wrong in a specific direction — concordance understated, and understated most in exactly the repetitive regions where callers actually differ. The tools here compare properly.
Section 10.6 is the reason to read their output stratified rather than as a single figure.
Its vcfeval command compares two callsets by haplotype rather than by position, so two spellings of the same indel are scored as agreement instead of as one false positive and one false negative. Position-wise comparison systematically understates concordance in exactly the regions where callers differ most, which makes this the difference between a benchmark and a number. (Cleary et al. 2015)
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| hap.py | dormant | C++py |
The wrapper most benchmarking papers and the PrecisionFDA challenges report through, which makes it the format other people’s numbers arrive in. It stratifies results by genomic region, which is the part worth having: one overall F-score hides that most of the disagreement lives in a small fraction of the genome. |
11.4 What is not here
Somatic calling is adjacent to this page and is not surveyed on it. Tumour samples break the assumption every tool above is built on — that a site has one of a small number of genotypes at known frequencies — because a tumor is a mixture of clones at unknown proportions with unknown ploidy, contaminated by normal tissue. That is a different inference problem rather than a harder setting for the same one, and it deserves its own treatment.
Structural variants are absent for a related reason, covered when that part is written: the evidence for a rearrangement is not the pile of bases at a position but the pattern of how reads fail to align, which is a different question again.
11.5 Language coverage
Two of the callers here are Java, two are C++, and the two built around learned models are Python wrapping compiled code — which is the pattern to expect wherever a tool has a model to load and preprocessing to do fast.
The chips report the same thing they do elsewhere: what you can drive each tool from. As with alignment, the answer is that you drive them from the command line and read the VCF back in, and the language coverage that matters for this part of the work is the coverage on the file-format side rather than the caller side.