A Cautionary Note on the Use of Genotype Callers in Phylogenomics

Abstract

Next-generation-sequencing genotype callers are commonly used in studies to call variants from newly sequenced species. However, due to the current availability of genomic resources, it is still common practice to use only one reference genome for a given genus, or even one reference for an entire clade of a higher taxon. The problem with traditional genotype callers, such as the one from GATK, is that they are optimized for variant calling at the population level. However, when these callers are used at the phylogenetic level, the consequences for downstream analyses can be substantial. Here, we performed simulations to compare the performance between the genotype callers of GATK and ATLAS, and present their differences at various phylogenetic scales. We show that the genotype caller of GATK substantially underestimates the number of variants at the phylogenetic level, but not at the population level. We also found that the accuracy of heterozygote calls declines with increasing distance to the reference genome. We quantified this decline and found that it is very sharp in GATK, while ATLAS maintains high accuracy even at moderately divergent species from the reference. We further suggest that efforts should be taken towards acquiring more reference genomes per species, before pursuing high-scale phylogenomic studies. [ATLAS; efficiency of SNP calling; GATK; heterozygote calling; next-generation sequencing; reference genome; variant calling.].

Figure 1

Base-cases example trees: a birth--death tree A (left), and a recent-burst tree B (right). For each tree we generated three arbitrary rescalings: Large, Medium, and Small, each representing various types of divergence found in phylogenetic studies. The chosen reference genome is shown with an arrow. The colors of the tip labels are chosen to represent the increasing distance to the reference (these colors will be used again in Figs. 4 and 5). The command lines and parameters used to generate these trees are described in the Supplementary Material, Section A.1 available on Dryad.

Figure 2

Summary of the main steps taken in this study. The order of steps is given by the numbers in each box. Initial sequence simulation and indexing of the reference genome involve haploid sequences (steps 1, 2, 4, and 6). Simulation of reads, mapping, recalibration, and genotype calling involve diploid sequences (steps 3, 5, 7, 8, and 9). Arrows indicate that some output of a previous step will be used as input for a next step. Programs used, and a summary of the command lines are indicated in italics. A complete description of all command lines is given in the Supplementary Material, Section A available on Dryad.

Figure 3

Number of called variants for the tips of tree A (first row), and tree B (second row). The tips on the x-axis are ordered according to their distance to the reference (the reference being always at the extreme left).

Figure 4

Accuracy of heterozygote calling for tree A at three different phylogenetic scales: Large (first column), Medium (middle column), and Small (last column). The phylogenetic distance from each tip of the phylogeny to the reference is plotted against the heterozygote-calling accuracy (first row), and against the called versus true heterozygotes (second row). The colors of each point or symbol correspond to the tip colors shown in Fig. 1.

Figure 5

Accuracy of heterozygote calling for tree B at three different phylogenetic scales: Large (first column), Medium (middle column), and Small (last column). The phylogenetic distance from each tip of the phylogeny to the reference is plotted against the heterozygote-calling accuracy (first row), and against the called versus true heterozygotes (second row). The colors of each point or symbol correspond to the tip colors shown in Fig. 1.