It is tempting to think of a variant caller as a program that finds SNPs. It is more accurate — and much more useful — to think of it as a program that answers a narrower question, one site at a time:

Given the pile of reads covering this position, what is the probability of the data under each possible genotype?

Everything else — the VCF, the QUAL column, the filtering you will do afterward — is bookkeeping built on top of that one calculation. If you understand the calculation, the rest of the pipeline stops being magic.

Notice what the question does not ask. It does not ask "is this a real SNP." It asks "how well does each candidate genotype explain the reads I observed." The difference matters: the caller is making a probabilistic statement about your data, not a claim about biological truth. A site can have a beautiful genotype likelihood and still be an artifact, because the reads themselves were placed in the wrong location by the aligner. A variant caller cannot fix a mapping error. It can only, at best, notice that something is off.

Every base in a FASTQ file carries a Phred quality score:

$$ Q = -10 \log_{10}(\epsilon) $$

where $$\epsilon$$ is the probability that the base call is wrong. So $$Q30$$ means $$\epsilon = 0.001$$, a 1-in-1000 chance of error. Inverting:

$$ \epsilon = 10^{-Q/10} $$

This is the atom of the whole system. The caller treats each base as an independent observation with a known error rate, and asks how surprising that observation is under each genotype.

Consider a diploid individual at one position. There are three possible genotypes if we restrict ourselves to two alleles, reference A and alternate G: A/A, A/G, G/G.

Now suppose we see 10 reads: 4 carry A, 6 carry G.

For each genotype we compute $$P(\text{reads} \mid \text{genotype})$$ — the genotype likelihood.

Under A/A (homozygous reference): every G we observed must be a sequencing error. If each G base has error probability $$\epsilon$$, seeing six of them is very unlikely. The likelihood is tiny.

Under G/G (homozygous alternate): every A must be an error. Also unlikely — four errors is a lot.

Under A/G (heterozygous): each read is drawn from one of the two chromosomes with probability $$1/2$$. Seeing a 4:6 split is entirely unremarkable. The likelihood is high.

Formally, for a diploid genotype $$G = {a_1, a_2}$$ and a set of $$n$$ independent read bases $$b_1 \ldots b_n$$:

$$ P(\text{reads} \mid G) = \prod_{i=1}^{n} \frac{1}{2}\left[P(b_i \mid a_1) + P(b_i \mid a_2)\right] $$

where $$P(b_i \mid a)$$ is $$1 - \epsilon_i$$ if the base matches allele $$a$$, and $$\epsilon_i/3$$ if it does not (dividing by 3 because an error could produce any of the other three bases). The $$\frac{1}{2}$$ and the sum encode the coin flip over which of the two homologous chromosomes the read came from.

This model — independent reads, Phred-scaled error rates, a sum over the alleles in the genotype — is the statistical core of bcftools mpileup and is described in Li (2011), the foundational paper for the pileup-based approach [doc9].

A naive approach would be "call a het if the alternate allele fraction is between 0.2 and 0.8." Compare the two:

The likelihood framework naturally handles depth and base quality together. Two reads and forty reads produce the same allele fraction but wildly different confidence — and the likelihood captures that, while a fraction cutoff does not. This is the single most important conceptual advantage, and it is the reason you should propagate likelihoods downstream rather than hard genotypes whenever your coverage is low.

A likelihood is not yet a genotype call. To get from $$P(\text{reads} \mid G)$$ to $$P(G \mid \text{reads})$$ — what we actually want — we need Bayes' theorem:

$$ P(G \mid \text{reads}) = \frac{P(\text{reads} \mid G) , P(G)}{\sum_{G'} P(\text{reads} \mid G') , P(G')} $$

The term $$P(G)$$ is the prior: what we believed about this genotype before looking at the reads. And here is where non-model work starts to diverge from the textbook.

By default the prior encodes "most sites are not variable." bcftools call uses a mutation-rate prior with a default of $$1.1 \times 10^{-3}$$ — a value derived from human heterozygosity [doc9]. That default is a statement about human genomes, and you are probably not working on humans.

This is not a trivial detail:

• Inbred lines, selfing plants, clonal fungi, inbred mapping populations have far less heterozygosity than the default assumes. The prior is too permissive and you will over-call heterozygotes from noise.
• Outbred trees, wild insect populations, highly diverse crop landraces can have heterozygosity an order of magnitude higher than human. The prior is too conservative and you will systematically under-call real heterozygous sites, especially at low depth.

bcftools call -P FLOAT exposes this directly: lower it for inbred material, raise it for diverse outbred populations [doc9]. Very few people in your field adjust it. Some of them should.

When you call several samples jointly, the caller can use the allele frequency across samples to inform the prior at each site. A site where 8 of 10 individuals clearly carry the alternate allele makes it much more plausible that the ambiguous 9th individual is a real heterozygote rather than noise.

This is why joint calling beats calling each sample separately and merging VCFs afterward, and it is most valuable exactly where you need it most — low coverage. It also introduces a subtle trap: your genotype calls now depend on which samples you included in the cohort. Add or remove individuals and the calls can change. This is correct Bayesian behavior, not a bug, but it means you must document your cohort composition for reproducibility.

Open any VCF and you will see a FORMAT field like GT:AD:DP:GQ:PL. Three of those are the statistical output described above.

CHROM  POS   REF ALT QUAL  ... FORMAT        SAMPLE_1
scaf3  4412  A   G   184   ... GT:AD:DP:GQ:PL  0/1:12,14:26:99:412,0,358

PL — Phred-scaled genotype likelihoods. This is the most information-dense field in the file and the one students most often skip. The three numbers correspond to genotypes 0/0, 0/1, 1/1, and are normalized so the most likely genotype is 0:

$$ \text{PL}_G = -10 \log_{10} \frac{P(\text{reads} \mid G)}{P(\text{reads} \mid \hat{G})} $$

So 412,0,358 means: heterozygous is the best explanation; homozygous reference is $$10^{41.2}$$ times less likely; homozygous alternate is $$10^{35.8}$$ times less likely. A confident call.

Contrast with PL = 30,0,25. Still called het, but the alternatives are only a few hundred-fold less likely. At low coverage you will see a great deal of this.

GQ — genotype quality. Simply the second-smallest PL: the Phred-scaled confidence that the called genotype is correct rather than the runner-up. GQ 99 is capped and very confident; GQ 3 means the caller is essentially guessing.

QUAL — site quality. The Phred-scaled probability that there is no variant at this site. Note the distinction that trips everyone up: QUAL is about the site being variable, GQ is about this sample's genotype being right. You can have a rock-solid site (QUAL 2000) with a garbage genotype for one low-coverage individual (GQ 4).

AD — allelic depths. Raw counts per allele. Always request this (-a FORMAT/AD in bcftools mpileup) [doc9] — you need it for filtering, for allelic-balance checks, and for spotting paralog problems.

For any downstream analysis that can accept them — site frequency spectra, $$F_{ST}$$, admixture, GWAS — pass genotype likelihoods rather than hard genotype calls. Hard-calling at 5× throws away the caller's own honest admission of uncertainty and biases you toward homozygotes, because at low depth you frequently sample only one of the two chromosomes. Tools built around the likelihood framework (ANGSD and its relatives) exist precisely for this reason and are widely used in non-model population genomics, where 3–10× coverage across many individuals is often a better use of a budget than 30× on a few.

Before comparing specific programs, understand the split, because it explains nearly every performance difference you will read about [doc9].

Position-based (pileup) callers — bcftools mpileup | bcftools call — trust the aligner. At each reference position they tally the reads, compute genotype likelihoods as described above, and call under a Bayesian model. Fast, transparent, no training data required, works at any ploidy [doc9].

Local-reassembly / haplotype callers — GATK HaplotypeCaller, DeepVariant, DRAGEN — discard the aligner's local placement in regions that look active. They assemble candidate haplotypes from the reads, realign every read against those haplotypes, and only then genotype [doc9]. The indel is represented once, on the assembled haplotype, rather than being smeared inconsistently across reads.

That distinction is why reassembly callers win on indels, clustered variants, and complex regions [doc9]. It is also why they are slower and, in pathological high-depth or repeat-rich regions, can abort or misbehave [doc9].

bcftools mitigates its structural disadvantage with BAQ (Base Alignment Quality), which downweights base qualities near likely misalignments so a shaky column does not produce a confident false SNP. BAQ is on by default and materially changes your calls [doc9]. It suppresses false positives; it does not reconstruct the true indel.

bcftools mpileup -Ou -f reference.fa \
    -q 20 -Q 20 \
    -a FORMAT/DP,FORMAT/AD,FORMAT/SP \
    sample1.bam sample2.bam sample3.bam \
  | bcftools call -mv -Oz -o variants.vcf.gz
bcftools index variants.vcf.gz

Strengths. Fast, few moving parts, no reference resources beyond the FASTA, and — critically for this audience — it handles arbitrary ploidy [doc9]. It is the default recommendation for non-model organisms, organelles, microbes, and quick exploratory scans [doc9].

There is a specific structural reason it behaves well on messy genomes: bcftools caps the variant quality at the mapping quality. A read in a repeat gets a low MAPQ because it could have come from several places, and bcftools propagates that ambiguity into the variant score rather than ignoring it [doc7]. Reassembly callers assume the read was placed in roughly the right region and only needs local adjustment — an assumption that fails badly when the read belongs to a different copy of a transposon entirely [doc7].

Weaknesses. Materially worse on indels in homopolymers and short tandem repeats, and in low-mappability regions [doc9]. Not designed for somatic/low-VAF detection or for cohorts beyond roughly a hundred samples [doc9]. Produces no machine-learning quality score, so all filtering is manual [doc9].

Flags that change your results [doc9]:

The reference implementation of the reassembly approach, and the pipeline everything else is benchmarked against [doc9]. Its GVCF workflow (HaplotypeCaller -ERC GVCF → GenomicsDBImport → GenotypeGVCFs) is genuinely excellent engineering for scalable cohort genotyping, and the ability to add samples later without recalling everything is a real practical advantage for a growing project.

The problem for you is filtering, not calling.

GATK's recommended filtering method is VQSR, which trains a machine-learning model to separate true from false variants using curated truth sets — HapMap, 1000 Genomes, dbSNP. Those resources do not exist for the vast majority of non-human species [doc7]. This is not a minor inconvenience; it removes the entire filtering apparatus GATK was designed around.

The fallback is hard filtering with fixed thresholds (QD < 2, FS > 60, MQ < 40, MQRankSum < -12.5, ReadPosRankSum < -8) [doc7]. A simulated insect-population benchmark found that GATK's annotation scores largely failed to separate true positives from false positives, making principled threshold selection very difficult [doc6]. Applying the standard human-derived thresholds caused severe damage: recovery rates ranged from 81% down to 18%, depending on how diverged the reference genome was from the study population [doc6].

That last point deserves emphasis, because it describes the situation many of you are in. The study varied $$d$$, the divergence between the reference genome and the sampled population, from 0 to 2%. As $$d$$ increased, hard-filtering performance collapsed [doc6]. If you are mapping to a reference assembled from a different individual, population, subspecies, or congener — which is routine outside human genetics — this is your scenario.

The mechanism: at $$d = 0$$, filtering on MQ < 40 removed 80% of false positives while discarding only 5% of true positives. At $$d = 0.02$$, the same filter discarded 79% of true positives [doc6]. The threshold did not change; the biology it was tuned for did.

A useful diagnostic from that work: plot the distribution of QD scores. A bimodal distribution signals that a large fraction of your variants are fixed differences between your population and the reference, rather than segregating polymorphism — and that QD-based filtering will savage your recall [doc6][doc7]. A unimodal distribution means hard filtering is comparatively safe [doc7].

DeepVariant reframes the problem entirely. Rather than computing a likelihood from an explicit error model, it renders the read pileup at each candidate site as a multi-channel image tensor and classifies it with a convolutional neural network (InceptionV3) into three classes: homozygous reference, heterozygous, homozygous alternate [doc5].

The network learns error profiles and local haplotype structure directly from training data rather than having them hand-specified. This works remarkably well: on human WGS it reaches F1 ≈ 0.9961 for SNPs and 0.9959 for indels [doc5], and it has won multiple PrecisionFDA Truth Challenge categories [doc1]. It leads open-source tools on indels and difficult regions, and the accuracy gain is largest at low coverage [doc1][doc9].

Two hard constraints for this audience, straight from the documentation:

1. Diploid only. DeepVariant supports organisms where ploidy is two, because the supported genotype classes are exactly hom-ref, het, and hom-alt [doc1]. Polyploid plants and fungi are simply out of scope. No flag fixes this.

2. The shipped models are trained on human data only [doc1]. The developers explicitly maintain a guide on non-human calling pitfalls [doc1].

How much does the human-training issue actually hurt? Less than you might fear, but it varies. DeepVariant is documented as generalizing reasonably off its training sample [doc9], and it has been applied at scale to non-human data — the Rice 3K dataset was reanalyzed with it [doc3]. The developers' own work on mosquito pedigrees found that species-specific retrained models improved calling, assessed via rates of Mendelian violation [doc3]. Retraining is possible and documented [doc4], but it requires a truth set or a pedigree and real computational investment — realistically a project, not a lab exercise.

A practical note: choosing the right --model_type (WGS, WES, PACBIO, ONT_R104, HYBRID) matters, and using the wrong one degrades accuracy [doc9]. Cohort merging uses GLnexus, not GenotypeGVCFs [doc9].

A defensible default for most of you: diploid, Illumina, reference from your own or a close species, population-genomic goals → bcftools mpileup + call, filtered on QUAL and depth, with repeat regions excluded.

The benchmark evidence supports this. In the simulated insect population, bcftools produced 54–521× fewer false positives than GATK HaplotypeCaller in raw callsets, with comparable recovery rates, regardless of whether reads were mapped with bowtie2 or bwa-mem [doc6]. After filtering, bcftools achieved the highest F-scores in most conditions [doc6]. The authors' recommendation was explicit: bcftools mpileup should be the first choice for non-human studies [doc7].

Two caveats before you treat that as settled. First, it is a single simulation study of one insect genome, and simulated populations differ systematically from real ones — the authors say so themselves [doc7]. Second, it found one condition where GATK wins: at low coverage (5×), GATK with VQSR outperformed bcftools on F-score [doc6]. If you have shallow data and a way to build a training set, that is worth knowing.

Switch to DeepVariant when indels matter, when your organism is reliably diploid, and when you can afford to validate. Switch to GATK when you need the GVCF cohort workflow for a large and growing sample set, or when you have a reference closely matched to your population.

Nothing in this reading will cost you more than getting ploidy wrong, and nothing is easier to get wrong outside vertebrates.

The genotype model has to be able to represent your organism's genotypes. If it cannot, the caller does not warn you — it produces plausible-looking garbage [doc9].

• Calling a diploid as haploid halves your heterozygote sensitivity. Every het becomes a coin flip between the two homozygotes.
• Calling a haploid as diploid manufactures false heterozygous calls from every sequencing error and every paralogous mismapping [doc9].

Cases that come up constantly in this room:

Fungi. Many are haploid for most of their life cycle. --ploidy 1 [doc9]. Dikaryons are a genuinely hard case that neither simple setting handles honestly.

Polyploid plants. Autotetraploids have genotype classes AAAA, AAAB, AABB, ABBB, BBBB — a dosage series, not three states. Setting --ploidy 4 preserves that dosage rather than collapsing everything to "het" [doc9]. Allopolyploids are harder still: homeologous subgenomes look exactly like heterozygous sites to a caller, and the standard fix is subgenome-aware mapping, not a ploidy flag. Also note that DeepVariant is unavailable to you here [doc1].

Organelles. Chloroplast and mitochondrial genomes are effectively haploid, but true heteroplasmy is continuous in allele fraction — not 0, 0.5, or 1. No diploid or haploid genotype model can express it; real heteroplasmy work needs a somatic-style caller [doc9].

Pooled sequencing. Set --ploidy N to the number of chromosomes in the pool so allele counts are preserved [doc9].

# Haploid fungal isolate
bcftools mpileup -f ref.fa in.bam | bcftools call -m --ploidy 1 -Oz -o out.vcf.gz

# Autotetraploid
bcftools mpileup -f ref.fa in.bam | bcftools call -m --ploidy 4 -Oz -o out.vcf.gz

A raw VCF from any caller is an intermediate file [doc9].

bcftools norm -f reference.fa -m -any variants.vcf.gz -Oz -o norm.vcf.gz

Two VCFs can encode the identical biological haplotype with different records — an indel placed at either end of a homopolymer, or an MNP versus two split SNVs. Normalization left-aligns indels and splits multiallelic records so identical variants get identical representations. Do this before any comparison, annotation, or merge [doc9]. Skipping it is the most common cause of "my two pipelines disagree" when in fact they agree perfectly.

Since bcftools produces no VQSR or deep-learning score, filtering is on explicit annotations [doc9]. A reasonable starting point:

bcftools filter -e 'QUAL < 30 || FORMAT/DP < 5 || FORMAT/DP > 100' \
    -s LOWQUAL -Oz -o filtered.vcf.gz norm.vcf.gz

Note the upper depth bound. Sites with anomalously high coverage are collapsed repeats — two or more genome copies stacked into one reference locus — and they generate false heterozygotes with beautiful allele balance. This is one of the most reliable artifact signatures in non-model genomes, where assemblies are more fragmented and repeat-collapsed than the human reference. Set the ceiling from your own coverage distribution, not from a number in a paper.

The benchmark evidence suggests bcftools filtering is unusually forgiving: dropping variants with QUAL < 20 reduced false positives by 35–42% while costing under 1% of recovery [doc6]. That is a very good trade, and it reflects the fact that bcftools' QUAL is a genuinely calibrated probability rather than a heuristic score.

Across both callers in the insect benchmark, 79–85% of all false positives fell in repeats [doc6]. GATK was disproportionately affected, with false positives heavily concentrated in LTR and LINE transposable elements — a 514-fold excess over bcftools in TEs overall, and over 1300-fold in simple repeats and satellites [doc6].

The recommendation from that study is blunt: exclude variants in repeats before downstream analysis [doc7]. If you have a RepeatMasker annotation, use it:

bcftools view -T ^repeats.bed filtered.vcf.gz -Oz -o final.vcf.gz

If you do not have one, running RepeatMasker or building a quick repeat library is usually a better investment than tuning filter thresholds.

bcftools stats final.vcf.gz > stats.txt

Check the transition/transversion ratio. Ti/Tv is elevated above the null expectation of 0.5 in most organisms because of CpG deamination and the structure of the genetic code — around 2.0–2.1 genome-wide in humans, but species-specific, so do not import the human number as a target. What is diagnostic is a Ti/Tv near 0.5: that means your callset is dominated by random error, because sequencing errors have no transition bias.

Also check the heterozygosity per sample. An individual with double everyone else's het rate is usually contaminated or has a sample-swap problem, not interesting biology.

If you set out to compare two callers on your own data — a good exercise — there are three ways to get it wrong [doc9]:

1. Not normalizing first. Un-normalized records mismatch on representation alone and you will report differences that do not exist.
2. Using bcftools isec as your benchmark. A record-level set intersection overcounts both false positives and false negatives. Proper benchmarking is haplotype-aware (hap.py + vcfeval, within a confident-region BED) [doc9].
3. Reporting one genome-wide F1. The headline number hides everything that matters. On easy SNPs, every modern caller exceeds F1 0.999 [doc9]. The differences live in indels, repeats, and low-mappability regions. Stratify and report SNPs and indels separately.

And the deeper problem, specific to your situation: you have no truth set. Human benchmarking rests on GIAB, and nothing equivalent exists for green ash or your fungal isolate. Concordance between two pipelines is not accuracy — as Heng Li has pointed out, if two pipelines share an artifact that a third lacks, the third looks worse while actually being better [doc7]. Obtaining true variant positions genome-wide from empirical data is essentially impossible [doc7], which is exactly why the insect benchmark used simulation.

Practical substitutes when you lack a truth set:

• Mendelian violation rates, if you have a pedigree or a cross. This is what the DeepVariant team used for mosquitoes [doc3].
• Replicate concordance — sequence one individual twice and measure disagreement.
• Overlap with an independent genotyping platform, if a SNP array exists for your system.
• Ti/Tv and the site frequency spectrum as sanity checks — weak evidence individually, but a badly distorted SFS is a real warning sign.

1. A genotype likelihood is $$P(\text{reads} \mid \text{genotype})$$, built from Phred base qualities. It handles depth and quality jointly in a way that allele-fraction cutoffs cannot.
2. Calls combine likelihood with a prior, and the default prior is human-derived [doc9]. Know whether yours should be higher or lower.
3. PL is the honest output; GT is a lossy summary of it. Pass likelihoods downstream when coverage is low.
4. Pileup callers trust the aligner; reassembly callers re-derive the local alignment. That single difference explains most benchmark results [doc9].
5. For non-model diploids, bcftools is a strong and defensible default — fewer false positives, easier filtering, any ploidy [doc6][doc7][doc9]. DeepVariant leads on indels but is diploid-only and human-trained [doc1]. GATK's cohort workflow is excellent but its filtering strategy assumes resources you do not have [doc7].
6. Ploidy must match your biology or the output is silently wrong [doc9].
7. Repeats are where false positives live — 79–85% of them [doc6]. Excluding repeats will do more for your callset than any amount of threshold tuning.
8. You have no truth set. Design your validation around pedigrees, replicates, or independent platforms, and stay skeptical of your own callset.

The core methodological sources behind this reading:

• Li H. (2011) A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics 27(21):2987–2993 — the genotype-likelihood model underlying mpileup [doc9].
• Danecek P. et al. (2021) Twelve years of SAMtools and BCFtools. GigaScience 10(2):giab008 [doc9].
• DePristo M.A. et al. (2011) A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics 43(5):491–498 [doc9].
• Poplin R. et al. (2018) A universal SNP and small-indel variant caller using deep neural networks. Nature Biotechnology 36:983–987 [doc1][doc9].
• Lefouili M. & Nam K. (2022) The evaluation of Bcftools mpileup and GATK HaplotypeCaller for variant calling in non-human species. Scientific Reports 12:11331 — the central non-model benchmark cited throughout [doc6][doc7].
• Krusche P. et al. (2019) Best practices for benchmarking germline small-variant calls in human genomes. Nature Biotechnology 37:555–560 [doc9].
• Yun T. et al. (2021) Accurate, scalable cohort variant calls using DeepVariant and GLnexus. Bioinformatics [doc1].
• Google Genomics team, Improved non-human variant calling using species-specific DeepVariant models (2018) — mosquito pedigree analysis [doc3].