Interpreting Tajima's D
Tajima's D is easy to compute and easy to over-read. This page covers what
pixy's D does and does not tell you, what pixy counts when computing
it, and how to use it on real data.
The short version
D is sound for comparing windows, populations, or groups computed the same way on the same data. That is what it is for.
D is not a p-value. Its absolute value is not a test of neutrality: under neutrality its expectation is not 0 and its standard deviation is not 1.
pixycounts mutations, not segregating sites. On biallelic data this is identical to the classic estimator; at multiallelic sites it departs by design from tools that count sites.Multiallelic sites in real data are enriched for genotyping artifacts, and θW and D weight them more heavily than π does. Filter accordingly.
Why D is not a p-value
Tajima's D is often taught as "0 under neutrality, negative after a sweep,
positive under balancing selection," with |D| > 2 as a rule of thumb for
significance. That framing relies on assumptions that essentially never hold for
whole-genome resequencing data, and it is not specific to pixy — it applies
to every implementation of the statistic.
The expectation is not 0. Tajima's D is a ratio whose denominator is itself a random variable, so E[D] ≠ 0 even under a strictly neutral, infinite-sites coalescent. Tajima (1989) used a beta approximation for significance precisely because of this. Separately, real sequence data violates the infinite-sites assumption: when a site is hit by mutation more than once, both π and θW undercount, and D settles at a small negative finite-sites floor rather than at 0. That floor deepens with missingness and depends on θ.
The standard deviation is not 1. The variance term (Tajima 1989) is derived for a single non-recombining locus. Real windows contain many recombining genealogies, which average out and shrink the spread of D substantially.
Note
Versions 2.1.0 through 2.2.1 additionally inflated the denominator whenever missing data made the per-site observed-allele count ragged, shrinking |*D*| toward 0 — badly enough to matter (a true D of −1.84 was reported as −1.38 at 20% missing genotypes). This was a regression of the fix for issue #160 and is corrected; see Changelog. Data with no missing genotypes was never affected. If you have Tajima's D from an affected version computed on data with missing genotypes, recompute it.
In pixy's polyploid validation (θ ≈ 0.038, 10 kb windows, JC69 mutation) the
measured values were:
Regime |
mean D |
sd D |
|---|---|---|
No recombination, no recurrent mutation |
≈ −0.11 |
≈ 0.87 |
No recombination, recurrent mutation |
≈ −0.08 to −0.14 |
≈ 0.83 |
Recombination at r = μ (realistic) |
finite-sites floor |
≈ 0.19 |
These are not universal constants — they depend on θ, window size, recombination rate, sample size, and missingness. They are reported here to make one point: the quantity you would need to divide by to turn D into a z-score is not 1, is not knowable from the formula, and varies several-fold across regimes.
What this means in practice. Use D as a relative measure. Comparing D
across windows within a dataset, or between populations processed identically, is
exactly what the statistic supports and what pixy is built for. Reading a
single window's D = −0.4 as evidence of a sweep is not supported. If you need
calibrated significance, you need a null simulated to match your data's θ,
ploidy, missingness, and per-window recombination rate — an analysis pixy
does not perform for you.
What pixy counts
Watterson's estimator is derived from the number of mutations on the genealogy: E[η] = anθ follows from Poisson mutation along the coalescent tree. The familiar "number of segregating sites" is a shorthand for that quantity, exact only under strict infinite sites, where every mutation creates a new site and the two counts coincide.
A site carrying three or four alleles is, by definition, a place where that
shorthand breaks. pixy therefore counts a site with k observed alleles as
k − 1 mutations — Tajima's s*, the minimum (parsimony) number of mutations
per site (Tajima 1996) — in both θW and the D variance term.
Consequences worth knowing:
On biallelic data nothing changes. Every biallelic site has k − 1 = 1, so
pixyreduces exactly to the classic estimator and agrees withscikit-allel.On multiallelic data ``pixy`` deliberately differs from
scikit-allel,vcftools, and other site-counting implementations. This is a design choice, not a discrepancy to reconcile. Report it in your methods.A site fixed for the alternate allele is not segregating. Every sample being homozygous
1/1means one observed allele and zero mutations, even though the site is a variant relative to the reference.pixycounts it as 0. (Versions through 2.2.1 incorrectly counted such sites as segregating; see Changelog.)θ:sub:`W` carries a small downward bias under recurrent mutation.
s*is a parsimony minimum: homoplasy and back-mutation are invisible to it. In the validation regime above this left θW ≈ 4% below 4Neμ, flat across ploidy. No model-free estimator can remove this; doing so requires committing to a specific mutation model.
Multiallelic sites on real data
pixy's multiallelic validation simulates a finite-sites Jukes–Cantor process,
in which every third or fourth allele is a genuine recurrent mutation. Real
multiallelic calls are not so clean. Collapsed paralogs, mismapping, and
indel-adjacent miscalls all manufacture spurious extra alleles, and they do so
preferentially in the repetitive regions where mapping is hardest.
This matters more for θW and D than for π. Because a k-allele site contributes k − 1 mutations, a spurious third allele contributes twice what a spurious second allele would, whereas any single site's contribution to π is bounded. The validation establishes that the estimators are correct given the genotypes; it does not establish that a given set of multiallelic calls has earned that trust.
How much will this affect me?
The size of the multiallelic contribution scales with θ and with ploidy. From
pixy's θ sweep, here is how much diversity biallelic-only π misses — a
useful proxy for how much multiallelic sites matter in your data:
θ |
diploid |
octoploid |
|---|---|---|
0.005 |
−0.8% |
−0.6% |
0.01 |
−1.5% |
−2.0% |
0.025 |
−3.8% |
−5.4% |
0.05 |
−7.1% |
−10.3% |
0.10 |
−12.5% |
−17.6% |
At human-like diversity (θ ≈ 0.001) multiallelic sites are a rounding error; enable the flag and move on. At high diversity, in polyploids, or with large samples, they matter enough to be worth the care described below.
A practical checklist
Run both modes and diff them. Run
pixywith and without--include_multiallelic_snps. The difference is the multiallelic contribution. If it is far larger than the table above predicts for your θ, your multiallelic calls are telling you about your pipeline, not your population.Check the allele spectrum. Genuine recurrent mutation produces mostly 3-allele sites and rarely 4-allele ones. A pile-up of 4-allele sites, or multiallelic sites clustered in repeats or beside indels, is a mapping signature rather than a biological one.
Filter multiallelic sites at least as hard as biallelic ones. Depth caps, mapping quality, and repeat masking all apply. The k − 1 weighting means a spurious triallelic site costs you double.
Compare the two modes across windows, not their absolute values.
Reporting in a methods section
Because pixy's multiallelic θW and D differ by design from
site-counting implementations, state what you used:
Watterson's θ and Tajima's D were computed with
pixyvX, which counts mutations (Tajima'ss*, the minimum mutation count Σ(k − 1)) rather than segregating sites. On biallelic data this is identical to the standard estimators; at multiallelic sites it departs by design from implementations that count sites (e.g.scikit-allel,vcftools). θ estimates carry an irreducible ~4% downward bias under recurrent mutation.
References
Tajima, F. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123: 585–595.
Tajima, F. (1996). The amount of DNA polymorphism maintained in a finite population when the neutral mutation rate varies among sites. Genetics 143: 1457–1465.
Watterson, G.A. (1975). On the number of segregating sites in genetical models without recombination. Theoretical Population Biology 7: 256–276.
Bailey, N., Stevison, L. & Samuk, K. (2025). Correcting for bias in estimates of θW and Tajima's D from missing data in next-generation sequencing. Molecular Ecology Resources e14104.
Roychoudhury, A. & Wakeley, J. (2010). Sufficiency of the number of segregating sites in the limit under finite-sites mutation. Theoretical Population Biology 78: 118–122.