Alnus

Moorjani et al. (2011) (link) first observed that pairwise LD measurements across a panel of SNPs can be combined to enable accurate inference of the age of admixture, n. The crux of their approach was to harness the fact that the ALD between two sites x and y scales as e^−nd multiplied by the product of allele frequency differences δ(x)δ(y) in the mixing populations. While the allele frequency differences δ(⋅) are usually not directly computable, they can often be approximated. Thus, Moorjani et al. (2011) (link) formulated a method, ROLLOFF, that dates admixture by fitting an exponential decay e^−nd to correlation coefficients between LD measurements and surrogates for δ(x)δ(y). Note that Moorjani et al. (2011) (link) define z(x, y) as a sample correlation coefficient, analogous to the classical LD measure r, as opposed to the sample covariance (Equation 1) we use here; we find the latter more mathematically convenient.
We build upon these previous results by deriving exact formulas for weighted sums of ALD under a variety of weighting schemes that serve as useful surrogates for δ(x)δ(y) in practice. These calculations will allow us to interpret the magnitude of weighted ALD to obtain additional information about admixture parameters. Additionally, the theoretical development will generally elucidate the behavior of weighted ALD and its applicability in various phylogenetic scenarios.
Following Moorjani et al. (2011) (link), we partition all pairs of SNPs (x, y) into bins of roughly constant genetic distance, where ε is a discretization parameter inducing a discretization on d. Given a choice of weights w(⋅), one per SNP, we define the weighted LD at distance d as
Assume first that our weights are the true allele frequency differences in the mixing populations, i.e., w(x) = δ(x) for all x. Applying Equation 3, where F₂(A, B) is the expected squared allele frequency difference for a randomly drifting neutral allele, as defined in Reich et al. (2009) (link) and Patterson et al. (2012) (link). Thus, a(d) has the form of an exponential decay as a function of d, with time constant n giving the date of admixture.
In practice, we must compute an empirical estimator of a(d) from a finite number of sampled genotypes. Say we have a set of m diploid admixed samples from population C indexed by i = 1, …, m, and denote their genotypes at sites x and y by x_i, y_i ε {0, 1, 2}. Also assume we have some finite number of reference individuals from A and B with empirical mean allele frequencies ${\widehat{p}}_A(\cdot )$ and ${\widehat{p}}_B(\cdot )$ . Then our estimator is where $\widehat{\text{cov}(X,Y)}=\frac{1}{m-1}{\displaystyle \sum _{i=1}^m\left(x_i-\stackrel{-}{x}\right)\left(y_i-\stackrel{-}{y}\right)}$ is the usual unbiased sample covariance, so the expectation over the choice of samples satisfies $E\left[\widehat{a}\left(d\right)\right]=a\left(d\right)$ (assuming no background LD, so the ALD in population C is independent of the drift processes producing the weights).
The weighted sum is a natural quantity to use for detecting ALD decay and is common to our weighted LD statistic $\widehat{a}\left(d\right)$ and previous formulations of ROLLOFF. Indeed, for SNP pairs (x, y) at a fixed distance d, we can think of Equation 3 as providing a simple linear regression model between LD measurements z(x, y) and allele frequency divergence products δ(x)δ(y). In practice, the linear relation is made noisy by random sampling, as noted above, but the regression coefficient 2αβe^−nd can be inferred by combining measurements from many SNP pairs (x, y). In fact, the weighted sum in the numerator of Equation 5 is precisely the numerator of the least-squares estimator of the regression coefficient, which is the formulation of ROLLOFF given in Moorjani et al. (2012, Note S1). Note that measurements of z(x, y) cannot be combined directly without a weighting scheme, as the sign of the LD can be either positive or negative; additionally, the weights tend to preserve signal from ALD while depleting contributions from other forms of LD.
Up to scaling, our ALDER formulation is roughly equivalent to the regression coefficient formulation of ROLLOFF (Moorjani et al. 2012 , Note S1). In contrast, the original ROLLOFF statistic (Patterson et al. 2012 (link)) computed a correlation coefficient between z(x, y) and w(x)w(y) over . However, the normalization term in the denominator of the correlation coefficient can exhibit an unwanted d-dependence that biases the inferred admixture date if the admixed population has undergone a strong bottleneck (Moorjani et al. 2012 , Note S1) or in the case of recent admixture and large sample sizes. Beyond correcting the date bias, the $\widehat{a}\left(d\right)$ curve that ALDER computes has the advantage of a simple form for its amplitude in terms of meaningful quantities, providing us additional leverage on admixture parameters. Additionally, we will show that $\widehat{a}\left(d\right)$ can be computed efficiently via a new fast Fourier transform-based algorithm.

As noted above, our ALDER formulation of weighted LD hones the original two-reference admixture dating technique of ROLLOFF (Moorjani et al. 2011 (link)), correcting a possible bias (Moorjani et al. 2012 , Note S1), and the one-reference technique (Pickrell et al. 2012 (link)), improving statistical power. Pickrell et al. (2012) (link) also observed that weighted LD can be used to estimate ancestral mixing fractions. We further develop this application now.
The main idea is to treat our expressions for the amplitude of the weighted LD curve as equations that can be solved for the ancestry fractions α and β = 1 − α. First consider two-reference weighted LD. Given samples from an admixed population C and reference populations A′ and B′, we compute the curve $\widehat{a}\left(d\right)$ and fit it as an exponential decay plus affine term: $\widehat{a}\left(d\right)\approx \widehat{M}{e}^{-nd}+\widehat{K}$ . Let ${\widehat{a}}_0:=\widehat{M}+\widehat{K}/2$ denote the amplitude of the curve. Then Equation 10 gives us a quadratic equation that we can solve to obtain an estimate $\widehat{\alpha }$ of the mixture fraction α, $2\widehat{\alpha }\left(1-\widehat{\alpha }\right)F_2{\left(A^{″},B^{″}\right)}^2={\widehat{a}}_0,$ assuming we can estimate F₂(A″, B″)². Typically the branch-point populations A″ and B″ are unavailable, but their F₂ distance can be computed by means of an admixture tree (Lipson et al. 2012 ; Patterson et al. 2012 (link); Pickrell and Pritchard 2012 (link)). A caveat of this approach is that α and 1 − α produce the same amplitude and cannot be distinguished by this method alone; additionally, the inversion problem is ill-conditioned near α = 0.5, where the derivative of the quadratic vanishes.
The situation is more complicated when using the admixed population as one reference. First, the amplitude relation from Equation 11 gives a quartic equation in $\widehat{\alpha }$ : $2\widehat{\alpha }\left(1-\widehat{\alpha }\right){\left[\widehat{\alpha }{F}_2\left(A,R^{″}\right)-\left(1-\widehat{\alpha }\right)F_2\left(B,R^{″}\right)\right]}^2={\widehat{a}}_0.$ Second, the F₂ distances involved are in general not possible to calculate by solving allele frequency moment equations (Lipson et al. 2012 ; Patterson et al. 2012 (link)). In the special case that one of the true mixing populations is available as a reference, however—i.e., R′ = A—Pickrell et al. (2012) (link) demonstrated that mixture fractions can be estimated much more easily. From Equation 7, the expected amplitude of the curve is 2αβ³F₂(A, B)². On the other hand, assuming no drift in C since the admixture, allele frequencies in C are given by weighted averages of allele frequencies in A and B with weights α and β; thus, the squared allele frequency differences from A to B and C satisfy $F_2\left(A,C\right)={\beta }^2{F}_2\left(A,B\right),$ and F₂(A, C) is estimable directly from the sample data. Combining these relations, we can obtain our estimate $\widehat{\alpha }$ by solving the equation $2\widehat{\alpha }/\left(1-\widehat{\alpha }\right)={\widehat{a}}_0/F_2{\left(A,C\right)}^2.$ In practice, the true mixing population A is not available for sampling, but a closely related population A′ may be. In this case, the value of $\widehat{\alpha }$ given by Equation 12 with A′ in place of A is a lower bound on the true mixture fraction α (see Appendix A for theoretical development and Results for simulations exploring the tightness of the bound). This bounding technique is the most compelling of the above mixture fraction inference approaches, as prior methods cannot perform such inference with only one reference population. In contrast, when more references are available, moment-based admixture tree-fitting methods, for example, readily estimate mixture fractions (Lipson et al. 2012 ; Patterson et al. 2012 (link); Pickrell and Pritchard 2012 (link)). In such cases we believe that existing methods are more robust than LD-based inference, which suffers from the degeneracy of solutions noted above; however, the weighted LD approach can provide confirmation based on a different genetic mechanism.

We fit discretized weighted LD curves $\widehat{a}\left(d\right)$ as $\widehat{M}{e}^{-nd}+\widehat{K}$ from Equation 9, using least-squares to find best-fit parameters. This procedure is similar to ROLLOFF, but ALDER makes two important technical advances that significantly improve the robustness of the fitting. First, ALDER directly estimates the affine term K that arises from the presence of subpopulations with differing ancestry percentages by using interchromosome SNP pairs that are effectively at infinite genetic distance (Appendix A). The algorithmic advances we implement in ALDER enable efficient computation of the average weighted LD over all pairs of SNPs on different chromosomes, giving $\widehat{K}$ and, importantly, eliminating one parameter from the exponential fitting. In practice, we have observed that ROLLOFF fits are sometimes sensitive to the maximum inter-SNP distance d to which the weighted LD curve is computed and fit; ALDER eliminates this sensitivity.
Second, because background LD is present in real populations at short genetic distances and confounds the ALD signal (interfering with parameter estimates or producing spurious signal entirely), it is important to fit weighted LD curves starting only at a distance beyond which background LD is negligible. ROLLOFF used a fixed threshold of d > 0.5 cM, but some populations have longer-range background LD (e.g., from bottlenecks), and moreover, if a reference population is closely related to the test population, it can produce a spurious weighted LD signal due to recent shared demography. ALDER therefore estimates the extent to which the test population shares correlated LD with the reference(s) and fits only the weighted LD curve beyond this minimum distance as in our test for admixture (Appendix B).
We estimate standard errors on parameter estimates by performing a jackknife over the autosomes used in the analysis, leaving out each in turn. Note that the weighted LD measurements from individual pairs of SNPs that go into the computed curve $\widehat{a}\left(d\right)$ are not independent of each other; however, the contributions of different chromosomes can reasonably be assumed to be independent.