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

Mason Liang; Mikhail Shishkin; Anastasia Mikhailova; Vladimir Shchur; Rasmus Nielsen

doi:10.1371/journal.pgen.1010281

Abstract

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.

Author summary

We establish a theoretical framework to model 3-locus admixture linkage disequilibrium of an admixed population taking into account the effects of genetic drift, migration and recombination. The theory is used to develop a method for estimating the times of multiple admixtures events. We demonstrate the accuracy of the method on simulated data and we apply it to previously published data from Mexican and Colombian populations to explore the complex history of American populations in the post-Colombian period.

Citation: Liang M, Shishkin M, Mikhailova A, Shchur V, Nielsen R (2022) Estimating the timing of multiple admixture events using 3-locus linkage disequilibrium. PLoS Genet 18(7): e1010281. https://doi.org/10.1371/journal.pgen.1010281

Editor: Garrett Hellenthal, University College London, UNITED KINGDOM

Received: November 9, 2021; Accepted: June 2, 2022; Published: July 15, 2022

Copyright: © 2022 Liang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: Software is openly available in the repository https://github.com/Genomics-HSE/LaNeta. Data sharing is not applicable as this study only analyses previously published and openly available data. The analysis uses genomic data from Yoruba, Colombian and Mexican populations. Whole-genome sequence data are available from the third phase of 1000 Genome Project (http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/). All filters are included in LaNeta in the file prepararion.sh (https://github.com/Genomics-HSE/LaNeta/blob/main/utilites/preparation.sh).

Funding: MS and VS worked on this paper within the framework of the HSE University Basic Research Program (hse.ru). AM was supported by the grant RFBR 20-29-01028 (rfbr.ru). RN was supported by NIH grant R01GM138634 (nih.gov). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

There are many methods for inferring the presence of admixture, e.g. methods using simple summary statistics detecting deviations from phylogenetic symmetry [1–3] and methods 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–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–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–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 CHROMOPAINTER [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].

Description of the method

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, (1) 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′ to the right of y. The three-loci LD, as defined by [29] and [30], is given by (2)

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₃ 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 , 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 (3)

By taking expectation, we obtain

Consider an arbitrary number N of sites S₁, …, S_N. We assume that for any i. Then we have (4)

Hence, we conclude that (5)

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: so that cov(A_x, A_y, A_z) = (1, −1, −1, −1, 2)v₃. There is one entry in v₃ 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 v_3(T) = (1 − m_T, (1 − m_T)², (1 − m_T)², (1 − m_T)², (1 − m_T)³)′. The probabilities for subsequent generations can be found by left-multiplying 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 (6)

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

This is because (1, −1, −1, −1, 2) is a left eigenvector of both L and U, with corresponding eigenvalues (1 − 1/2N)(1 − 2/2N) and exp(−d − d′). 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₂ generations before the present, so that , and all other m_i’s are 0. Making the substitution T₁ = T − T₂, the right hand side of Eq 6 becomes (7)

The corresponding expression for the two-point covariance function is given by (8) 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′. 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₁ × S_n and n₂ × S_n, where n₁ and n₂ 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′ 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 where and similarly for and .

Algorithm

Directly computing over the set d, d′ ∈ {0, 1, …, P}² would be cubic in the number of segregating sites. However, by using the fast Fourier Transform (FFT) technique introduced in ALDER [26], we can approximate with an algorithm whose time complexity is instead linear in the number of segregating sites.

First, rearrange to get and define sequences b_i[d] and c[d] by binning the data and then doubling the length by padding with P zeros,

We can approximate |S[d, d′]| and the n sums in the numerator of 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′)-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² 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₁ and M₂, a similar equation holds (9)

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₁, T₂ 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 non-linear 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). The exception is when the total admixture proportion M₂ + M₁(1 − M₂) is close to 0.5. As the total admixture proportion increases above 0.5, the contours for Eq 2 flip from being concave down to concave up. This transition can be seen by comparing the upper left side of Fig 2 to its lower right. At this threshold, the contours of the estimated weighted LD depend on the actual admixture fractions of the samples, which may differ from the expectation as a result of genetic drift. This mismatch between theory and simulations is most evident in Fig 2, for m₁ = 0.1, m₂ = 0.4 and m₁ = 0.2, m₂ = 0.4.

Download:

Fig 1. Predicted weighted LD surfaces from simulations and theory for varying admixture times.

The heat maps are from simulations and the contours are plotted from Eq 7. The two admixture probabilities were fixed at m₁ = m₂ = .2 and the the times of the two admixture pulses, T₁ and T₂, were varied. Each square covers the range 0.5 cM < d, d′ < 20 cM. When time of the more recent pulse is greater than half of that of the more ancient pulse, i.e. 2T₁ > T₁ + T₂, the contours of the resulting weighted LD surface are straight, making it difficult to distinguish from the weighted LD surface produced by a one-pulse admixture scenario.

https://doi.org/10.1371/journal.pgen.1010281.g001

Download:

Fig 2. Predicted weighted LD surfaces from simulations and theory for varying admixture proportions.

The heat maps are from simulations and the contours are plotted from Eq 7. The two admixture times were fixed at 2 and 12 generations ago (T₁ = 10 and T₂ = 2) while the admixture probabilities were varied. Each square covers the range 0.5 cM < d, d′ < 20 cM. As the total admixture proportion m₂ + m₁(1 − m₂) increases above 0.5, the contours change to reflecting that the majority contribution of the genetic material now originates from the other population. Weighted LD surfaces for m₁ > 0.5 or m₂ > 0.5 are not shown, but are qualitatively similar to the surfaces on the lower and rightmost sides.

https://doi.org/10.1371/journal.pgen.1010281.g002

Download:

Fig 3. Weighted LD surfaces produced by constant admixture.

The heat maps are from simulations and the contours from analytical results for a model in which continuous admixture started 10, 20, or 40 generations ago and stopped 5 generations before the present. Each square covers the range 0.5 cM < d, d′ < 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

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₁ > T₂). 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 log-scaled LD curves. Notice (Fig 4) that the slopes are somewhat similar for admixture models with identical values of the most recent admixture events (T₂). 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).

Download:

Fig 4. Two-locus weighted LD with two admixture events and varying pulse times.

Corresponding ALDER curves for two-pulse admixture with varying pulse times. Morgans on x-axis and log ALDER scores on y-axis. Red lines are T₂ = 2, Green lines T₂ = 5, and blue lines are T₂ = 10.

https://doi.org/10.1371/journal.pgen.1010281.g004

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₁ = 0 in our models). There is less variance in this estimate, as it can be explained by a single unknown parameter (T₂) compared to three free parameters estimated in our method (T₁, T₂ and admixture proportion M₁).

Download:

Fig 5. Accuracy of estimates of T₁ (A) and T₂ (B), and ALDER estimates of admixture time (C) as a function of other parameters.

Twelve admixture scenarios, T₁ ∈ {0, 5, 10, 20} and T₂ ∈ {2, 5, 10}, were simulated 100 times each. The admixture probabilities were fixed at M₁ = 0.3 and M₂ = 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₁ = 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₂ > T₁. ALDER estimates a single admixture time which corresponds to T₁ = 0.

https://doi.org/10.1371/journal.pgen.1010281.g005

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.

Download:

Fig 6. Effect of admixture proportion misspecification on the estimated values of T₁ and T₂.

Admixture proportion misspecification has a strong effect on the estimates of time T₁ between pulses of admixture. Estimates of the time T₂ of the most recent admixture pulse remain stable.

https://doi.org/10.1371/journal.pgen.1010281.g006

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 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].

Download:

Fig 7. Weighted LD surface for Mexican samples with Yoruba as the first source population reference.

The model with the best fit is two pulses from the non-Yoruba source population at T₁ + T₂ = 13.2 ± 1.01 and T₂ = 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 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.

Download:

Fig 8. Weighted LD surface for Colombian samples with Yoruba as the first source population reference.

The two-pulse model that fits best is two pulses of non-Yoruba admixture at T₁ + T₂ = 14.5 ± 0.74 and T₂ = 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

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 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.

Acknowledgments

This research was supported in part through computational resources of HPC facilities at HSE University [39].

References

1. Reich D, Thangaraj K, Patterson N, Price AL, Singh L. Reconstructing Indian population history. Nature. 2009;461(7263):489–494. pmid:19779445
- View Article
- PubMed/NCBI
- Google Scholar
2. Patterson N, Moorjani P, Luo Y, Mallick S, Rohland N, Zhan Y, et al. Ancient admixture in human history. Genetics. 2012;192(3):1065–93. pmid:22960212
- View Article
- PubMed/NCBI
- Google Scholar
3. Durand EY, Patterson N, Reich D, Slatkin M. Testing for ancient admixture between closely related populations. Mol Biol Evol. 2011;28(8):2239–52. pmid:21325092
- View Article
- PubMed/NCBI
- Google Scholar
4. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155(2):945–959. pmid:10835412
- View Article
- PubMed/NCBI
- Google Scholar
5. Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Research. 2009;19(9):1655–1664. pmid:19648217
- View Article
- PubMed/NCBI
- Google Scholar
6. Maples BK, Gravel S, Kenny EE, Bustamante CD. RFMix: A Discriminative Modeling Approach for Rapid and Robust Local-Ancestry Inference. The American Journal of Human Genetics. 2013;93(2):278–288. pmid:23910464
- View Article
- PubMed/NCBI
- Google Scholar
7. Pool JE, Nielsen R. Inference of historical changes in migration rate from the lengths of migrant tracts. Genetics. 2009;181(2):711–719. pmid:19087958
- View Article
- PubMed/NCBI
- Google Scholar
8. Gravel S. Population genetics models of local ancestry. Genetics. 2012;191(2):607–619. pmid:22491189
- View Article
- PubMed/NCBI
- Google Scholar
9. Liang M, Nielsen R. The Lengths of Admixture Tracts. Genetics. 2014; p. genetics–114. pmid:24770332
- View Article
- PubMed/NCBI
- Google Scholar
10. Corbett-Detig R, Nielsen R. A Hidden Markov Model Approach for Simultaneously Estimating Local Ancestry and Admixture Time Using Next Generation Sequence Data in Samples of Arbitrary Ploidy. PLOS Genetics. 2017;13(1):1–40. pmid:28045893
- View Article
- PubMed/NCBI
- Google Scholar
11. Svedberg J, Shchur V, Reinman S, Corbett-Detig R. Inferring Adaptive Introgression Using Hidden Markov Models. Molecular Biology and Evolution. 2021;38. pmid:33502512
- View Article
- PubMed/NCBI
- Google Scholar
12. Ni X., Yuan K., Yang X., Feng Q., Guo W., Ma Z. & Xu S. Inference of multiple-wave admixtures by length distribution of ancestral tracks. Heredity. 121, 52–63 (2018) pmid:29358727
- View Article
- PubMed/NCBI
- Google Scholar
13. Fisher RA. The Theory of Inbreeding. Edinburgh, Scotland: Oliver and Boyd; 1949.
14. Stam P. The distribution of the fraction of the genome identical by descent in finite random mating populations. Genetics Research. 1980;35:131–155.
- View Article
- Google Scholar
15. Guo SW. Computation of identity-by-descent proportions shared by two siblings. American Journal of Human Genetics. 1994;54(6):1104. pmid:8198133
- View Article
- PubMed/NCBI
- Google Scholar
16. Bickeböller H, Thompson EA. Distribution of genome shared IBD by half-sibs: approximation by the Poisson clumping heuristic. Theoretical Population Biology. 1996;50(1):66–90. pmid:8776838
- View Article
- PubMed/NCBI
- Google Scholar
17. Bickeböller H, Thompson EA. The probability distribution of the amount of an individual’s genome surviving to the following generation. Genetics. 1996;143(2):1043–1049. pmid:8725249
- View Article
- PubMed/NCBI
- Google Scholar
18. Stefanov VT. Distribution of genome shared identical by descent by two individuals in grandparent-type relationship. Genetics. 2000;156(3):1403–1410. pmid:11063711
- View Article
- PubMed/NCBI
- Google Scholar
19. Ball F, Stefanov VT. Evaluation of identity-by-descent probabilities for half-sibs on continuous genome. Mathematical Biosciences. 2005;196(2):215–225. pmid:15979108
- View Article
- PubMed/NCBI
- Google Scholar
20. Cannings C. The identity by descent process along the chromosome. Human heredity. 2003;56(1-3):126–130. pmid:14614246
- View Article
- PubMed/NCBI
- Google Scholar
21. Dimitropoulou P, Cannings C. RECSIM and INDSTATS: probabilities of identity in general genealogies. Bioinformatics. 2003;19(6):790–791. pmid:12691999
- View Article
- PubMed/NCBI
- Google Scholar
22. Walters K, Cannings C. The probability density of the total IBD length over a single autosome in unilineal relationships. Theoretical Population Biology. 2005;68(1):55–63. pmid:15927222
- View Article
- PubMed/NCBI
- Google Scholar
23. Rodolphe F, Martin J, Della-Chiesa E. Theoretical description of chromosome architecture after multiple back-crossing. Theoretical Population Biology. 2008;73(2):289–299. pmid:18190938
- View Article
- PubMed/NCBI
- Google Scholar
24. Baird SJ, Barton NH, Etheridge AM. The distribution of surviving blocks of an ancestral genome. Theoretical Population Biology. 2003;64(4):451–471. pmid:14630482
- View Article
- PubMed/NCBI
- Google Scholar
25. Moorjani P, Patterson N, Hirschhorn JN, Keinan A, Hao L, Atzmon G, et al. The history of African gene flow into Southern Europeans, Levantines, and Jews. PLoS genetics. 2011;7(4):e1001373. pmid:21533020
- View Article
- PubMed/NCBI
- Google Scholar
26. Loh PR, Lipson M, Patterson N, Moorjani P, Pickrell JK, Reich D, et al. Inferring admixture histories of human populations using linkage disequilibrium. Genetics. 2013;193(4):1233–1254. pmid:23410830
- View Article
- PubMed/NCBI
- Google Scholar
27. Moreno-Mayar JV, Rasmussen S, Seguin-Orlando A, Rasmussen M, Liang M, Flåm ST, et al. Genome-wide Ancestry Patterns in Rapanui Suggest Pre-European Admixture with Native Americans. Current Biology. 2014;. pmid:25447991
- View Article
- PubMed/NCBI
- Google Scholar
28. Gravel S, Zakharia F, Moreno-Estrada A, Byrnes JK, Muzzio M, Rodriguez-Flores JL, et al. Reconstructing native American migrations from whole-genome and whole-exome data. PLoS genetics. 2013;9(12):e1004023. pmid:24385924
- View Article
- PubMed/NCBI
- Google Scholar
29. Bennett J. On the theory of random mating. Annals of Eugenics. 1952;17(1):311–317.
- View Article
- Google Scholar
30. Slatkin M. On treating the chromosome as the unit of selection. Genetics. 1972;72(1):157–168. pmid:4672513
- View Article
- PubMed/NCBI
- Google Scholar
31. Hellenthal G., Busby G., Band G., Wilson J., Capelli C., Falush D. & Myers Simon A Genetic Atlas of Human Admixture History. Science. 343, 747–751 (2014), https://www.science.org/doi/abs/10.1126/science.1243518 pmid:24531965
- View Article
- PubMed/NCBI
- Google Scholar
32. Lawson D., Hellenthal G., Myers S. & Falush D. Inference of population structure using dense haplotype data. PLoS Genetics. 8, e1002453 (2012) pmid:22291602
- View Article
- PubMed/NCBI
- Google Scholar
33. Liang M, Nielsen R. Understanding admixture fractions. bioRxiv. 2014; p. 008078.
- View Article
- Google Scholar
34. Verdu P, Rosenberg NA. A general mechanistic model for admixture histories of hybrid populations. Genetics. 2011;189(4):1413–1426. pmid:21968194
- View Article
- PubMed/NCBI
- Google Scholar
35. Pickrell JK, Pritchard JK. Inference of population splits and mixtures from genome-wide allele frequency data. PLoS genetics. 2012;8(11):e1002967. pmid:23166502
- View Article
- PubMed/NCBI
- Google Scholar
36. Kelleher J, Etheridge AM, McVean G. Efficient coalescent simulation and genealogical analysis for large sample sizes. PLoS computational biology. 2016;12(5):e1004842. pmid:27145223
- View Article
- PubMed/NCBI
- Google Scholar
37. Pedersen BS, Quinlan AR. cyvcf2: fast, flexible variant analysis with Python. Bioinformatics. 2017;. pmid:28165109
- View Article
- PubMed/NCBI
- Google Scholar
38. Consortium GP, et al. A global reference for human genetic variation. Nature. 2015;526(7571):68.
- View Article
- Google Scholar
39. Kostenetskiy P, Chulkevich R, Kozyrev V. HPC Resources of the Higher School of Economics. In: Journal of Physics: Conference Series. vol. 1740. IOP Publishing; 2021. p. 012050.

[ref1] 1. Reich D, Thangaraj K, Patterson N, Price AL, Singh L. Reconstructing Indian population history. Nature. 2009;461(7263):489–494. pmid:19779445
View Article
PubMed/NCBI
Google Scholar

[2] View Article

[3] PubMed/NCBI

[4] Google Scholar

[ref2] 2. Patterson N, Moorjani P, Luo Y, Mallick S, Rohland N, Zhan Y, et al. Ancient admixture in human history. Genetics. 2012;192(3):1065–93. pmid:22960212
View Article
PubMed/NCBI
Google Scholar

[6] View Article

[7] PubMed/NCBI

[8] Google Scholar

[ref3] 3. Durand EY, Patterson N, Reich D, Slatkin M. Testing for ancient admixture between closely related populations. Mol Biol Evol. 2011;28(8):2239–52. pmid:21325092
View Article
PubMed/NCBI
Google Scholar

[10] View Article

[11] PubMed/NCBI

[12] Google Scholar

[ref4] 4. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155(2):945–959. pmid:10835412
View Article
PubMed/NCBI
Google Scholar

[14] View Article

[15] PubMed/NCBI

[16] Google Scholar

[ref5] 5. Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Research. 2009;19(9):1655–1664. pmid:19648217
View Article
PubMed/NCBI
Google Scholar

[18] View Article

[19] PubMed/NCBI

[20] Google Scholar

[ref6] 6. Maples BK, Gravel S, Kenny EE, Bustamante CD. RFMix: A Discriminative Modeling Approach for Rapid and Robust Local-Ancestry Inference. The American Journal of Human Genetics. 2013;93(2):278–288. pmid:23910464
View Article
PubMed/NCBI
Google Scholar

[22] View Article

[23] PubMed/NCBI

[24] Google Scholar

[ref7] 7. Pool JE, Nielsen R. Inference of historical changes in migration rate from the lengths of migrant tracts. Genetics. 2009;181(2):711–719. pmid:19087958
View Article
PubMed/NCBI
Google Scholar

[26] View Article

[27] PubMed/NCBI

[28] Google Scholar

[ref8] 8. Gravel S. Population genetics models of local ancestry. Genetics. 2012;191(2):607–619. pmid:22491189
View Article
PubMed/NCBI
Google Scholar

[30] View Article

[31] PubMed/NCBI

[32] Google Scholar

[ref9] 9. Liang M, Nielsen R. The Lengths of Admixture Tracts. Genetics. 2014; p. genetics–114. pmid:24770332
View Article
PubMed/NCBI
Google Scholar

[34] View Article

[35] PubMed/NCBI

[36] Google Scholar

[ref10] 10. Corbett-Detig R, Nielsen R. A Hidden Markov Model Approach for Simultaneously Estimating Local Ancestry and Admixture Time Using Next Generation Sequence Data in Samples of Arbitrary Ploidy. PLOS Genetics. 2017;13(1):1–40. pmid:28045893
View Article
PubMed/NCBI
Google Scholar

[38] View Article

[39] PubMed/NCBI

[40] Google Scholar

[ref11] 11. Svedberg J, Shchur V, Reinman S, Corbett-Detig R. Inferring Adaptive Introgression Using Hidden Markov Models. Molecular Biology and Evolution. 2021;38. pmid:33502512
View Article
PubMed/NCBI
Google Scholar

[42] View Article

[43] PubMed/NCBI

[44] Google Scholar

[ref12] 12. Ni X., Yuan K., Yang X., Feng Q., Guo W., Ma Z. & Xu S. Inference of multiple-wave admixtures by length distribution of ancestral tracks. Heredity. 121, 52–63 (2018) pmid:29358727
View Article
PubMed/NCBI
Google Scholar

[46] View Article

[47] PubMed/NCBI

[48] Google Scholar

[ref13] 13. Fisher RA. The Theory of Inbreeding. Edinburgh, Scotland: Oliver and Boyd; 1949.

[ref14] 14. Stam P. The distribution of the fraction of the genome identical by descent in finite random mating populations. Genetics Research. 1980;35:131–155.
View Article
Google Scholar

[51] View Article

[52] Google Scholar

[ref15] 15. Guo SW. Computation of identity-by-descent proportions shared by two siblings. American Journal of Human Genetics. 1994;54(6):1104. pmid:8198133
View Article
PubMed/NCBI
Google Scholar

[54] View Article

[55] PubMed/NCBI

[56] Google Scholar

[ref16] 16. Bickeböller H, Thompson EA. Distribution of genome shared IBD by half-sibs: approximation by the Poisson clumping heuristic. Theoretical Population Biology. 1996;50(1):66–90. pmid:8776838
View Article
PubMed/NCBI
Google Scholar

[58] View Article

[59] PubMed/NCBI

[60] Google Scholar

[ref17] 17. Bickeböller H, Thompson EA. The probability distribution of the amount of an individual’s genome surviving to the following generation. Genetics. 1996;143(2):1043–1049. pmid:8725249
View Article
PubMed/NCBI
Google Scholar

[62] View Article

[63] PubMed/NCBI

[64] Google Scholar

[ref18] 18. Stefanov VT. Distribution of genome shared identical by descent by two individuals in grandparent-type relationship. Genetics. 2000;156(3):1403–1410. pmid:11063711
View Article
PubMed/NCBI
Google Scholar

[66] View Article

[67] PubMed/NCBI

[68] Google Scholar

[ref19] 19. Ball F, Stefanov VT. Evaluation of identity-by-descent probabilities for half-sibs on continuous genome. Mathematical Biosciences. 2005;196(2):215–225. pmid:15979108
View Article
PubMed/NCBI
Google Scholar

[70] View Article

[71] PubMed/NCBI

[72] Google Scholar

[ref20] 20. Cannings C. The identity by descent process along the chromosome. Human heredity. 2003;56(1-3):126–130. pmid:14614246
View Article
PubMed/NCBI
Google Scholar

[74] View Article

[75] PubMed/NCBI

[76] Google Scholar

[ref21] 21. Dimitropoulou P, Cannings C. RECSIM and INDSTATS: probabilities of identity in general genealogies. Bioinformatics. 2003;19(6):790–791. pmid:12691999
View Article
PubMed/NCBI
Google Scholar

[78] View Article

[79] PubMed/NCBI

[80] Google Scholar

[ref22] 22. Walters K, Cannings C. The probability density of the total IBD length over a single autosome in unilineal relationships. Theoretical Population Biology. 2005;68(1):55–63. pmid:15927222
View Article
PubMed/NCBI
Google Scholar

[82] View Article

[83] PubMed/NCBI

[84] Google Scholar

[ref23] 23. Rodolphe F, Martin J, Della-Chiesa E. Theoretical description of chromosome architecture after multiple back-crossing. Theoretical Population Biology. 2008;73(2):289–299. pmid:18190938
View Article
PubMed/NCBI
Google Scholar

[86] View Article

[87] PubMed/NCBI

[88] Google Scholar

[ref24] 24. Baird SJ, Barton NH, Etheridge AM. The distribution of surviving blocks of an ancestral genome. Theoretical Population Biology. 2003;64(4):451–471. pmid:14630482
View Article
PubMed/NCBI
Google Scholar

[90] View Article

[91] PubMed/NCBI

[92] Google Scholar

[ref25] 25. Moorjani P, Patterson N, Hirschhorn JN, Keinan A, Hao L, Atzmon G, et al. The history of African gene flow into Southern Europeans, Levantines, and Jews. PLoS genetics. 2011;7(4):e1001373. pmid:21533020
View Article
PubMed/NCBI
Google Scholar

[94] View Article

[95] PubMed/NCBI

[96] Google Scholar

[ref26] 26. Loh PR, Lipson M, Patterson N, Moorjani P, Pickrell JK, Reich D, et al. Inferring admixture histories of human populations using linkage disequilibrium. Genetics. 2013;193(4):1233–1254. pmid:23410830
View Article
PubMed/NCBI
Google Scholar

[98] View Article

[99] PubMed/NCBI

[100] Google Scholar

[ref27] 27. Moreno-Mayar JV, Rasmussen S, Seguin-Orlando A, Rasmussen M, Liang M, Flåm ST, et al. Genome-wide Ancestry Patterns in Rapanui Suggest Pre-European Admixture with Native Americans. Current Biology. 2014;. pmid:25447991
View Article
PubMed/NCBI
Google Scholar

[102] View Article

[103] PubMed/NCBI

[104] Google Scholar

[ref28] 28. Gravel S, Zakharia F, Moreno-Estrada A, Byrnes JK, Muzzio M, Rodriguez-Flores JL, et al. Reconstructing native American migrations from whole-genome and whole-exome data. PLoS genetics. 2013;9(12):e1004023. pmid:24385924
View Article
PubMed/NCBI
Google Scholar

[106] View Article

[107] PubMed/NCBI

[108] Google Scholar

[ref29] 29. Bennett J. On the theory of random mating. Annals of Eugenics. 1952;17(1):311–317.
View Article
Google Scholar

[110] View Article

[111] Google Scholar

[ref30] 30. Slatkin M. On treating the chromosome as the unit of selection. Genetics. 1972;72(1):157–168. pmid:4672513
View Article
PubMed/NCBI
Google Scholar

[113] View Article

[114] PubMed/NCBI

[115] Google Scholar

[ref31] 31. Hellenthal G., Busby G., Band G., Wilson J., Capelli C., Falush D. & Myers Simon A Genetic Atlas of Human Admixture History. Science. 343, 747–751 (2014), https://www.science.org/doi/abs/10.1126/science.1243518 pmid:24531965
View Article
PubMed/NCBI
Google Scholar

[117] View Article

[118] PubMed/NCBI

[119] Google Scholar

[ref32] 32. Lawson D., Hellenthal G., Myers S. & Falush D. Inference of population structure using dense haplotype data. PLoS Genetics. 8, e1002453 (2012) pmid:22291602
View Article
PubMed/NCBI
Google Scholar

[121] View Article

[122] PubMed/NCBI

[123] Google Scholar

[ref33] 33. Liang M, Nielsen R. Understanding admixture fractions. bioRxiv. 2014; p. 008078.
View Article
Google Scholar

[125] View Article

[126] Google Scholar

[ref34] 34. Verdu P, Rosenberg NA. A general mechanistic model for admixture histories of hybrid populations. Genetics. 2011;189(4):1413–1426. pmid:21968194
View Article
PubMed/NCBI
Google Scholar

[128] View Article

[129] PubMed/NCBI

[130] Google Scholar

[ref35] 35. Pickrell JK, Pritchard JK. Inference of population splits and mixtures from genome-wide allele frequency data. PLoS genetics. 2012;8(11):e1002967. pmid:23166502
View Article
PubMed/NCBI
Google Scholar

[132] View Article

[133] PubMed/NCBI

[134] Google Scholar

[ref36] 36. Kelleher J, Etheridge AM, McVean G. Efficient coalescent simulation and genealogical analysis for large sample sizes. PLoS computational biology. 2016;12(5):e1004842. pmid:27145223
View Article
PubMed/NCBI
Google Scholar

[136] View Article

[137] PubMed/NCBI

[138] Google Scholar

[ref37] 37. Pedersen BS, Quinlan AR. cyvcf2: fast, flexible variant analysis with Python. Bioinformatics. 2017;. pmid:28165109
View Article
PubMed/NCBI
Google Scholar

[140] View Article

[141] PubMed/NCBI

[142] Google Scholar

[ref38] 38. Consortium GP, et al. A global reference for human genetic variation. Nature. 2015;526(7571):68.
View Article
Google Scholar

[144] View Article

[145] Google Scholar

[ref39] 39. Kostenetskiy P, Chulkevich R, Kozyrev V. HPC Resources of the Higher School of Economics. In: Journal of Physics: Conference Series. vol. 1740. IOP Publishing; 2021. p. 012050.

Figures

Abstract

Author summary

Introduction

Description of the method

Model

Linkage disequilibrium and local ancestry.

Local ancestry covariance functions.

Weighted linkage disequilibrium.

Algorithm

Missing source population

Fitting the two-pulse model.

Verification and comparison

Patterns of 3-locus LD

Comparison to two-locus LD measures

Accuracy of parameter estimates

Applications

Discussion

Acknowledgments

References