(C) PLOS One
This story was originally published by PLOS One and is unaltered.
. . . . . . . . . .



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

['Mason Liang', 'Department Of Integrative Biology', 'University Of California', 'Berkeley', 'California', 'United States Of America', 'Mikhail Shishkin', 'International Laboratory Of Statistical', 'Computational Genomics', 'Hse University']

Date: 2022-09

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 ith 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 ith 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 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 , 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 1 , …, 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 3 . 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 v 3(T) = (1 − m T , (1 − m T )2, (1 − m T )2, (1 − m T )2, (1 − m T )3)′. 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 2 generations before the present, so that , and all other m i ’s are 0. Making the substitution T 1 = T − T 2 , 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 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′ 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}2 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 jth 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(P2 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 (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 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 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.

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]. PPT PowerPoint slide

PNG larger image

TIFF original image 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 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 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. PPT PowerPoint slide

PNG larger image

TIFF original image 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 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

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

[END]
---
[1] Url: https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1010281

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/

You are viewing proxied material from magical.fish. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.