Lewontin's D' for two biallelic loci from unphased genotypes

Formula: with \(p_A, p_B\) the major-allele frequencies, \(p_a = 1-p_A\), \(p_b = 1-p_B\) and \(p_{AB}\) the haplotype frequency, \(D = p_{AB} - p_A p_B\), \(D_{max} = \min(p_A p_b,\, p_a p_B)\), \(D_{min} = \max(-p_A p_B,\, -p_a p_b)\), and \(D' = D/D_{max}\) if \(D > 0\) else \(D/D_{min}\). Normalising by the attainable extreme is the whole content of the proposal: raw D is bounded by the allele frequencies and cannot be compared across loci, whereas D' reaches 1 exactly when one haplotype is absent. Note the sign convention, which is the one genetics::LD uses and is easy to misread: \(D_{min}\) is negative, so dividing a negative D by it gives a positive D'. estimate is therefore the normalised magnitude in [0, 1] and the direction is carried by the sign of D; sources that put D' on [-1, 1] mean sign(D) * estimate. Genotypes are unphased, so only the double heterozygote is ambiguous; it is resolved by EM run for a fixed 500 iterations with no convergence test, so the answer is deterministic and identical in both arms.

Usage

Dprime(geno1, geno2)

Arguments

geno1, geno2

Genotypes coded as the allele count 0, 1, 2. Individuals miscoded or missing at either locus are dropped pairwise.

Value

List with estimate (D'), D, pAB, pA, pB, Dmax, Dmin, r, r2, n, method.

References

Lewontin (1964), Genetics 49:49-67. The Genetics PDF at PMC could not be retrieved from this host (the fetch returned a 1.8 kB error page), so the coded form was read from Warnes and Leisch's CRAN package genetics, R/LD.R, which gives Dmin, Dmax and estDp verbatim. genetics maximises the same likelihood with optimize(); a fixed-iteration EM is used here because a golden-section search is not reproducible across language arms.

Examples

V <- c(1, 2, 3, 4, 5, 6, 7, 8)
Dprime(V, V)
#> $estimate
#> [1] 1
#> 
#> $D
#> [1] 0.1875
#> 
#> $pAB
#> [1] 0.75
#> 
#> $pA
#> [1] 0.75
#> 
#> $pB
#> [1] 0.75
#> 
#> $Dmax
#> [1] 0.1875
#> 
#> $Dmin
#> [1] -0.0625
#> 
#> $r
#> [1] 1
#> 
#> $r2
#> [1] 1
#> 
#> $n
#> [1] 2
#> 
#> $method
#> [1] "Lewontin D' two-locus disequilibrium"
#> 
#> attr(,"class")
#> [1] "morie_rich_result" "list"