Estimating the timing of multiple admixture events using 3-locus linkage disequilibrium

Estimating admixture histories is crucial for understanding the genetic diversity we see in present-day populations. Allele frequency or phylogeny-based methods are excellent for inferring the existence of admixture or its proportions. However, to estimate admixture times, spatial information from admixed chromosomes of local ancestry or the decay of admixture linkage disequilibrium (ALD) is used. One popular method, implemented in the programs ALDER and ROLLOFF, uses two-locus ALD to infer the time of a single admixture event, but is only able to estimate the time of the most recent admixture event based on this summary statistic. To address this limitation, we derive analytical expressions for the expected ALD in a three-locus system and provide a new statistical method based on these results that is able to resolve more complicated admixture histories. Using simulations, we evaluate the performance of this method on a range of different admixture histories. As an example, we apply the method to the Colombian and Mexican samples from the 1000 Genomes project. The implementation of our method is available at https://github.com/Genomics-HSE/LaNeta.


Introduction
There are many methods for inferring the presence of admixture, e.g. methods using simple summary statistics detecting deviations from phylogenetic symmetry [1][2][3]  estimating admixture proportions using programs such as Structure [4], Admixture [5] or RFmix [6]. There has also been substantial research on estimating admixture times. Some approaches are based on inferring admixture tract length distributions, such as [7][8][9][10][11][12]. Over time, recombination is expected to decrease the average lengths of admixture tracts. The length distribution of admixture tracts is therefore informative about the time since admixture. Much of the theory relating to tracts lengths is based on Fisher's famous theory of junctions [13] and subsequent work, such as [14][15][16][17][18][19][20][21][22][23]. For example, [24] first discussed the length distribution of tracts descended from a single ancestor. These results informed later analyses of admixture tract length distribution, such as references [7][8][9]. Gravel [8] also implemented the software program TRACTS, which estimates admixture histories by fitting the tract length distribution, obtained by local ancestry inference, to a exponential approximation.
Another approach, which we will follow in this paper, is based on the decay of admixture linkage disequilibrium (ALD). Linkage disequilibrium exists in any natural population due to mutation and genetic drift. However, in well-mixed and genetically isolated populations with recombination it usually decays quite rapidly at a genomic scale. For example, in many human populations linkage disequilibrium decays to approx. zero in less than 1 Mb. However, admixture tracts introduced into a population in an admixture event generates ALD over much longer distances, even if the amount of LD in the source populations is negligible. After a single admixture event, linkage disequilibrium in the admixed population will then gradually decrease in the subsequent generations as a result of recombination. It is, therefore, possible to make inferences about the admixture history of a population from the patterns of LD present in the population. This insight was first used in the program ROLLOFF [25] and was later extended by ALDER [26]. These two methods use the fact that if an admixed population takes in no additional migrants after the founding generation, the ALD present in the population is expected to decay approximately exponentially as a function of distance. The rate constant of this exponential decay is proportional to the age of the founding admixture pulse and can be used as an estimator. ROLLOFF and ALDER are well suited for inferring the time of the admixture event when the admixture history of the population can be approximated as a single pulse. However, in many realistic scenarios the admixture histories involve multiple pulses. Prominent examples in humans include Native American admixture in Rapa Nui [27] or admixed population groups in the Americas [28]. In these instances the expected decay of LD will become a mixture of exponentials. Existing dating method based on ALD can usually only infer the date of the most recent migration wave [25], or reject the hypothesis of a single pulse admixture [26].
ROLLOFF and ALDER use the information contained in pairs of sites by examining the two-locus linkage disequilibrium between them. Here we extend the theory underlying the methods in ROLLOFF and ADLER to three loci by considering three-locus LD. There are two ways of measuring the linkage between n loci. Bennett [29] defines n-locus linkage in a way that maintains a geometric decrease of LD each generation as a result of recombination, which is an important property of two-locus linkage disequilibrium. Slatkin [30] defines n-locus LD to be the n-way covariance, analogously to the property of two locus LD as the covariance in allele frequency between pairs of loci. For two and three loci, these two definitions coincide, but for four or more loci, they do not. Another method GLOBETROTTER [31] uses a copying model in a way similar to ROLLOFF. First, shared haplotype chunks are inferred with CHRO-MOPAINTER [32], then a mixture of exponentials is used to fit the coancestry curve.
In this paper, we will use Bennett and Slatkin's definition of three-locus LD to examine the decay of ALD for three sites as a function of the genetic distance between them. We derive an equation that describes the decay of three-locus LD under an admixture history with multiple waves of migration from two source populations. We then compare the results of coalescent simulations to this equation, and develop some guidelines for when admixture histories more complex than a single pulse can be resolved using ALD. Finally, we apply our method to the Colombian and Mexican samples in the 1000 Genomes data set, using the Yoruba samples as a reference. Fitting a two-pulse model to data, we estimate admixture histories for the two populations which are qualitatively consistent with the results reported in [28].

Model
We use a random union of gametes admixture model as described in [33], which is an extension of the mechanistic admixture model formulated by [34]. In this model, two or more source populations contribute migrants to form an admixed population consisting of 2N haploid individuals. Each generation in the admixed population is formed through the recombination of randomly selected individuals from the previous generation, with some individuals potentially replaced by migrants from the source populations. For simplicity, we consider a model with only two source populations. Furthermore, the first source population only contributes migrants in the founding generation, T. The second source population contributes migrants in the founding generation and possibly in one or more generations thereafter. In generation i, for i = T − 1, . . ., 0 (before the present), a fraction m i of the admixed population is replaced by individuals from the second source population.
Linkage disequilibrium and local ancestry. ROLLOFF and ALDER use the standard two-locus measure of LD between a SNP at positions x and another SNP at position y, which is a genetic distance d to the right, where H x and H y represent the haplotype or genotypes of an admixed chromosome at positions x and y. In the case of haplotype data, H i,x = 1 if the i th sample is carrying the derived allele at the SNP at position x, and is otherwise 0. Alternatively, for genotype data, H i,x take on values from {0, 1/2, 1} depending on the number of copies of the derived allele the i th sample is carrying at SNP position x. We consider an additional site at position z, which is located a further genetic distance d 0 to the right of y. The three-loci LD, as defined by [29] and [30], is given by The LD in an admixed population depends on the genetic differentiation between the source populations and and its admixture history. Let A x represent the local ancestry at position x, with A x = 0 if x is inherited from an ancestor in the second source population (the one which contributed in two admixture events), and A x = 1 if x is inherited from the first source population (the one which contributed in a single admixture event). We can compute D 3 in terms of the three-point covariance function of A x and so separate out the effects of allele frequencies and local ancestry. Consider the conditional expectation EðH x jA where g x is the allele frequency of locus x in the second source population and δ x = f x − g x is the difference of the allele frequencies of locus x in the two source populations. We now make the assumption that the allele frequencies in the source populations are known and fixed. Our goal is to prove that By taking expectation, we obtain Consider an arbitrary number N of sites S 1 , . . ., S N . We assume that EðH S i jA S 1 ; . . . ; Hence, we conclude that In particular, we obtain Eq 3 with N = 3. Local ancestry covariance functions. From the above section we see that we can describe the three-point admixture LD in terms of covariances of local ancestry in the three points. We now expand the covariance in Eq 2 into its component expectations to get Each one of these expectations on the right-hand side is the probability that one or more sites is inherited from an ancestor from the first source population. We organize these products of probabilities in a column vector: There is one entry in v 3 for each of the five ways in which the three markers at positions x, y, and z can arranged on one or more chromosomes. In the founding generation T, this column vector is given by The probabilities for subsequent generations can be found by leftmultiplying drift, recombination, and migration matrices: The matrices D i , L, and U account for the effects of migration, drift, and recombination, respectively. The migration matrix is a diagonal matrix given by Its entries are the probabilities that one, two, or three chromosomes in the admixed population will not be replaced by chromosomes from the second source population in generation i. The lower triangular drift matrix gives the standard Wright-Fisher drift transition probabilities between the states as a function of the population size 2N. Finally, the upper triangular recombination matrix is determined by the recombination rates between the three sites: The covariance function is then given by We can obtain an analogous equation for cov(A x , A y ), involving the migration, drift, and recombination matrices for two loci: In some cases, Eq 6 simplifies further. In a one-pulse migration model, in which the admixture proportion in the founding generation is m T = M and is there after 0, the D i 's become identity matrices, and we get the closed from expression is a left eigenvector of both L and U, with corresponding eigenvalues (1 − 1/2N)(1 − 2/2N) and exp(−d − d 0 ). Note that when M = 0, the covariance function will be identically 0. Another case is a two pulse model in which we ignore the effects of genetic drift. In this model, admixture only occurs T and T 2 generations before the present, so that m T ¼ M 1 ; m T 2 ¼ M 2 , and all other m i 's are 0. Making the substitution T 1 = T − T 2 , the right hand side of Eq 6 becomes The corresponding expression for the two-point covariance function is given by which is a mixture of two exponentials. Weighted linkage disequilibrium. As [26] noted, we cannot use the LD in the admixed population directly, because the allele frequency differences in the source populations can be of either sign. As in [26], we solve this problem by computing the product of the values of the three-point linkage disequilibrium coefficient with the product of the allele frequency differences. Using Eq 3 we obtain because the local ancestry in the admixed sample is independent of the allele frequencies in the admixed population. For inference purposes, we estimate this function by averaging over triples of SNPs which are separated by distances of approximately d and d 0 . The LD term is estimated from the admixed population, while the δ's are estimated from reference populations which are closely related to the two source populations. We notice that both this approach, as well as the previous approaches (e.g., [26]), do not take genetic drift in the source populations after the time of admixture into account, i.e. there is an assumption of both this method and previous methods that the allele frequencies in the ancestral source populations can be approximated well using the allele frequencies in the extant populations. We arrange the data from the admixed samples in an n × S n matrix H, where n is the number of admixed haplotypes/genotypes, and S n is the number of markers in the sample. Similarly, we arrange the data from the two source populations into two matrices, F and G, which are of size n 1 × S n and n 2 × S n , where n 1 and n 2 are the numbers of samples from each of the source populations. For ease of notation, we assume that the positions are given in units which make the unit interval equal to the desired bin width.
For a given d and d 0 the SNP triples we use in the estimator for the weighted LD are Let h x be empirical allele frequency in the admixed population. An estimator of the weighted three-point linkage disequilibrium coefficient is then and similarly ford y andd z .

Algorithm
We can approximate |S[d, d 0 ]| and the n sums in the numerator ofâ½d; d 0 � in terms of convolutions of these sequences: These convolutions can be efficiently computed with an FFT, since under a two-dimensional discrete Fourier transform from (d, d 0 )-space to (j, k)-space, where B i is the one-dimensional discrete Fourier transform of b and for j > 0, B i [−j] is the j th to last most element of B i . Summing over i and taking the inverse discrete Fourier transform, we can approximate the discrete Fourier transform of the numerator ofâ. We apply the same method to c to approximate the denominator ofâ. The time complexities for the binning and the FFT's are O(S n ) and O(P 2 log(P)). Of these two, the first term will dominate, because P, the number of bins, is much smaller than S n , the number of segregating sites.

Missing source population
When data from only one source population are available, it is still possible to estimate the weighted admixture linkage disequilibrium by estimating the difference in allele frequencies between the two source populations using the allele frequency differences between the available population and the admixed population [26,35], by way of the following formula where M is the admixture proportion. For two pulses of admixture with proportions M 1 and M 2 , a similar equation holds The allele frequencies in the missing source population can be estimated from this equation by solving for the relevant unknown term (either f x or g x ). This estimator might be noisy for rare variants, so sites with minor allele frequencies of less than 0.05 should be removed (this corresponds to standard filtering practices for real data).
When using only the admixed population itself as a reference population, the method described above will be biased if the same samples are used to estimate both the linkage disequilibrium coefficients and the weights (δ x , δ y , and δ z ). We cannot efficiently compute a polyache statistics like [26]. At the cost of some power, we instead adopt the approach of [35] and separate the admixed population into two equal-sized groups. We then use one group to estimate the weights, and the other group to estimate linkage disequilibrium coefficients, and vice versa. This gives two unbiased estimates for the numerator ofâ, which we then average.
Another challenge with real data, is that the method might be unstable when admixture proportions are not known. For two pulses of admixture, we have four independent parameters T 1 , T 2 and M1, M2. In order to simplify the problem, one can estimate the total ancestry fractions of the source populations in the admixed population using ADMIXTURE [5] with K = 2. Assume that M is the ancestry fraction of the source population which admixed two times. Then This equations are closely related to Eq 9. This allows to reduce the number of independent parameters to 3, which simplifies the optimization problem substantially.
Fitting the two-pulse model. We fit Eq 7 to the estimates of the weighted LD using nonlinear least squares, with two modifications. We added a proportionality constant to account for the expected square allele frequency difference between the source populations. We also subtracted out an affine term in the weighted LD which is due to population substructure [26]. We estimated this by computing the three-way covariance between triples of chromosomes. We use the jackknife to obtain confidence intervals for the resulting estimates by leaving out each chromosome in turn and refitting on the data for the remaining chromosomes.

Verification and comparison
We used the package msprime [36] to generate two source populations which diverged 4000 generations ago and a coalescent simulation to generate an admixed population from the two source populations according to two-pulse and constant admixture models. We sampled 50 diploid individuals from the admixed and two source populations, each consisting of 20 chromosomes of length 1 Morgan. The effective population size was 2N = 1000 for the admixed population and two source populations. Using a two pulse model, we varied the migration probabilities and timings for each pulse to examine the accuracy of Eq 7. We also simulated data for a model with a constant rate of admixture each generation, and compared this to the predictions made by Eq 6.
Our implementation uses Python package cyvcf2 [37] to read VCF files.

Patterns of 3-locus LD
We first evaluate the accuracy of the equations developed in this paper by comparing the analytical results to simulated data (Figs 1-3). We find there is a generally a close match between our equations and the simulated data under both the two-pulse admixture scenarios (Figs 1  and 2) and constant-admixture scenarios (Fig 3). When there is continuous admixture scenario, the shape of the weighted LD surface depends on both the duration and total amount of admixture. When the duration is short, the weighted LD surfaces are indistinguishable from the weighted LD surfaces produced by one pulse of migration. As the duration increases, the contours of the weighted LD surface become more curved. The contours are concave up when the total proportion is greater than 50% and concave down when it is less. When the total proportion is exactly 50%, the amplitude of the weighted LD surface is much smaller than the sampling error.
For two pulse models, the effects of the second pulse of migration only become evident when temporal spacing between the pulses is large enough (T 1 > T 2 ). Otherwise, the resulting weighted LD surface cannot be distinguished from the weighted LD surface produced by one pulse of admixture. As in the case of continuous admixture the concavity of the surface contours is determined by the total admixture proportion.

Comparison to two-locus LD measures
We compared the simulation results to the two-locus weighted LD calculated by ALDER (Fig  4). The information used in estimating Admixture times in ALDER is the slope of the logscaled LD curves. Notice (Fig 4) that the slopes are somewhat similar for admixture models with identical values of the most recent admixture events (T 2 ). Hence, when two admixture events have occurred, estimation of admixture times tend to get weighted towards the most recent event. Generally, it would be very difficult, based on the shape of the admixture LD decay curve to estimate parameters of a model with more than one admixture event. In contrast, there is a quite clear change in the pattern of three-locus LD as long as the time between the two admixture events is sufficiently large (Fig 1).

Accuracy of parameter estimates
We next evaluate the utility of the method for estimating admixture times. The qualitative similarities between one pulse and two pulse admixture scenarios seen in the previous simulations under some parameter settings will naturally affect the estimates. As shown in Fig 5, when the spacing between the two pulses is small relative to their age, the median of the estimates of the timing of the second pulse is close to the true value, but the interquartile range is large. Moreover, the best fit often lies on a boundary of the parameter space which is equivalent to a one pulse admixture model. When the spacing between the pulses is larger, the estimates for the timing of the older pulse become more precise. ALDER estimates a single admixture time (which corresponds to T 1 = 0 in our models). There is less variance in this estimate, as it can be explained by a single unknown parameter (T 2 ) compared to three free parameters estimated in our method (T 1 , T 2 and admixture proportion M 1 ).
We evaluated how admixture proportion mis-specification might affect the admixture time estimates. The results are summarised in Fig 6. The timing of the most recent admixture pulse is rather stable to variation in the admixture proportion, while the timing of the older pulse turns out to be quite sensitive to it.

Applications
To illustrate the utility of the method we computed weighted LD surfaces for Mexican and Colombian samples individuals in the third phase of the 1000 Genomes data set [38]. These samples were previously analyzed for similar purposes by [28]. Our datasets consisted of 64 individuals from Los Angeles and 94 individuals from Medellin, respectively. We used the 104 Yoruba samples as a reference population. We removed indels and SNVs and leave SNPs that only refers to autosomes. (All filters are included in utilites/preparation.sh in https://github. com/Genomics-HSE/LaNeta.) We computed the weighted LD on the genotypes to avoid effects of phasing errors.
For the Mexican samples, [28] found a small but consistent amount of African ancestry, which appeared in the population 15 generations ago, with continuing contributions from European and Native American populations since that date, but no African migration. We fitted a two-pulse model to the Mexican weighted LD surface (Fig 7) with Yoruba as the first source population and the other population being modeled as a missing source population, as Each square covers the range 0.5 cM < d, d 0 < 20 cM. We varied the time of the beginning of the admixture and the total admixture probability. The admixture probability for each generation was constant, and chosen so that the total admixture proportion was either 0.3 or 0.7. When the admixture is spread over 5 generations (the leftmost column), the resulting weighted LD surface is similar to a one-pulse weighted LD surface. For longer durations, the weighted LD surfaces are similar to those produced by two pulses of admixture.
https://doi.org/10.1371/journal.pgen.1010281.g003 previously described. The missing source population (non-Yoruban) represents an unknown admixed European and Native American population. This model was chosen to mimic the previous analysis of reference [28] which used a similar model to approximate continued gene-flow from Europeans and Native Americans. Using this set-up we estimated that the two pulses occurred 13.2 ± 1.01 and 7.9 ± 0.99 generations ago. Our results are roughly consistent with those of [28].
The weighted LD surface for the Colombia samples with Yoruba as the first source population is shown in Fig 8. From this, we estimated two pulses of non-Yoruba migration at 14.5 ± 0.74 and 3.7 ± 0.62 generations before the present. [28] inferred two pulses of admixture, corresponding to 13 and 5 generations ago. The weighted LD surface of the Colombian samples has contours which are strongly concave up, in contrast to those of the Mexican samples.

Discussion
The method presented here is an extension of previously published methods for using weighted two-locus LD to estimate admixture times. The new method uses more information in the data because it compares triples of SNPs instead of pairs. This gives the method the ability to infer admixture histories more complex than a one-pulse model. However, this comes at the price of greater estimation variances. ALDER and ROLLOFF make estimates from just tens of samples, while our method requires hundreds of samples. Part of this difference can be attributed to the fact that ALDER and ROLLOFF make inferences over a smaller class of models, but the main reason arises from the fact that the two-locus methods are estimating second  The admixture probabilities were fixed at M 1 = 0.3 and M 2 = 0.2. The colored bars give the medians of estimates for each of these twelve cases, the boxes delimit the interquartile range, and the whiskers extend out to 1.5 times the interquartile range. As the time between the two pulses of admixture increases, the error in the estimates decreases (for this reason we do not include T 1 = 0 accuracy estimate, in this case the results become unreasonable). Consistent with the simulations shown in Fig 1, there is limited power to estimate the time of the more ancient admixture pulse when T 2 > T 1 . ALDER estimates a single admixture time which corresponds to T 1 = 0.
https://doi.org/10.1371/journal.pgen.1010281.g005  The model with the best fit is two pulses from the non-Yoruba source population at T 1 + T 2 = 13.2 ± 1.01 and T 2 = 7.9 ± 0.99 generations ago. The weighted LD surface was estimated from real data, the level lines correspond to the best-fitting model inferred by LaNeta method.
https://doi.org/10.1371/journal.pgen.1010281.g007 The twopulse model that fits best is two pulses of non-Yoruba admixture at T 1 + T 2 = 14.5 ± 0.74 and T 2 = 3.7 ± 0.62 generations ago. The amplitude of this weighted LD surface is approximately ten times larger than that of the Mexican samples. This is a result of larger proportion of Yoruba ancestry in the Colombian samples. The weighted LD surface was estimated from real data, the level lines correspond to the best-fitting model inferred by LaNeta method. https://doi.org/10.1371/journal.pgen.1010281.g008

PLOS GENETICS
moments of the data, while we are estimating third moments. The variance of these estimates are both inversely proportional to the sample size, but the constants for estimating third moments are larger. As data becomes more readily available, this disadvantage should disappear.
We also note that the theory developed in this paper might be useful for other purposes than estimating admixture times. In particular, it can be used to test hypotheses regarding the spatial distribution of introgressed fragments in the genome, without relying on particular inferences of admixture tracts.