Modeling Continuous Admixture Using Admixture-Induced Linkage Disequilibrium

Recent migrations and inter-ethnic mating of long isolated populations have resulted in genetically admixed populations. To understand the complex population admixture process, which is critical to both evolutionary and medical studies, here we used admixture induced linkage disequilibrium (LD) to infer continuous admixture events, which is common for most existing admixed populations. Unlike previous studies, we expanded the typical continuous admixture model to a more general scenario with isolation after a certain duration of continuous gene flow. Based on the new models, we developed a method, CAMer, to infer the admixture history considering continuous and complex demographic process of gene flow between populations. We evaluated the performance of CAMer by computer simulation and further applied our method to real data analysis of a few well-known admixed populations.

In the present study, we first developed a weighted LD-based method to infer admixture with HI, GA, and continuous gene flow (CGF) 13 models (see Fig. 1). Both GA and CGF models assume that gene flow is a continuous process. Next, we extended the GA and CGF models to GA-I and CGF-I models, respectively (see Fig. 1), which models a scenario with a continuous gene flow duration followed by a period of isolation to present. We applied our method to a number of well-known admixed populations and provided information that would help better understanding the admixture history of these populations.

Material and Methods
Datasets. Data for simulation and empirical analysis were obtained from three public resources: Human Genome Diversity Panel (HGDP) 14 , the International HapMap Project phase III 15 and the 1000 Genomes Project (1KG) 16 . Source populations for simulations are the haplotypes from 113 Utah residents with Northern and Western European ancestries from the CEPH collection (CEU) and 113 Africans from Yoruba (YRI).
Inferring Admixture Histories by Using HI, GA, and CGF Models. The expectation of weighted LD under a two-way admixture model has been described in detail in another work 11 . Following the previous notation, the expectation of weighted LD statistic between two sites separated by a distance d (in Morgan) is as follows:  12 2 , δ 12 (x) is the allele frequency difference between populations 1 and 2 at site x, and S(d) is the set holding pairs of SNPs of distance d; a i (d), i = 0, 1, 2 are the weighted LD statistics of the admixed population (i = 0) and the source population i, (i = 1, 2), respectively; m i is the admixture proportion from the source populationi; c (l) is admixture indicator for the admixture event of l generations ago, and n is supposed to be the number of generations since the source populations first met. To eliminate the confounding effect due to background LD from the source populations, we used the quantity, z(d), defined as follows, to represent the admixture induced LD (ALD) 11 .
For different admixture models where admixture began ngenerations ago, z(d) varies in terms of the vector of coefficients of polynomial functions 12 : where the vector C model has length n using the HI, GA, CGF1, or CGF2 model; and n represents when the admixture occurred (HI) or began (GA and CGF) in terms of generations. For different models, the coefficient vectors have different patterns (see Fig. 2), which can be used to infer the best-fit model for a certain admixed population. In the CGF model, CGF1 represents the admixture where source population 1 is the recipient of the gene flow from population 2, whereas CGF2 indicates source population 2 as gene flow recipient from population 1. Inference of the admixture time assuming the true admixture history in one of these different models can be regarded as minimizing the objective function as follows: The optimization problem is therefore expressed as follows: where Z = (z(d 1 ), z(d 2 ), … , z(d I )) T is the observed ALD calculated from the single nucleotide polymorphism (SNP) data of both the parental populations and the admixed population, both Z and admixture proportion m i can be calculated by the algorithm iMAAPs 12 ; θ 0 is a real number used to correct the population substructure; θ 1 is a scalar that improves estimation robustness; 1 ∈ R I is a vector with each entry being 1; A is an I × J matrix with the ith row vector defined as Poly(d i ) T , i.e., A = (Poly(d 1 ), Poly(d 2 ), … , Poly(d I )) T , and J ≥ n is a pre-specified upper bound of n. Our definitions are consistent since we can let all entries after the n-th entry be 0 in C model . Next, we tried to estimate the parameters θ 0 , θ 1 , and C model , where C model has the information of the admixture model and the related admixture time n (in generations). In our analysis, the value of n is assumed to be a positive integer; therefore, our method is to go through all possible n values (with a reasonable upper limit J) to estimate n with the minimum value of the objective function. Given n, we used the ordinary least squares method to estimate (θ 0 , θ 1 ) such that the objective function was minimized. Using this approach, the value of n in relation to the minimal objective function value for each model was determined, which represents the time of admixture occurrence under each model. The method to conclude which models are the best is described in Identification of the best-fit model session.
Admixture Inference under HI, GA-I, and CGF-I Models. GA and CGF models assume that the admixture is strictly continuous from the beginning of admixture to present. This assumption seems too strong to be valid in empirical studies. Here, we extended GA model and CGF model to GA-I model and CGF-I model respectively, by considering continuous admixture followed by isolation. In this case, the admixture event lasts from G start generations ago to G end generations ago. Similar to the previous case, the coefficients of polynomial functions can be represented as a vector of length G start for each model, whose first G end − 1 entries are filled with zeros. Suppose the admixture lasted for n = G start − G end + 1 generations, then In this case, we can also try to find the parameters that minimize the objective function (eq. 2) under new models. By examining all possible pairs of (G end , G start ), it is possible to determine the global minimum of the objective function, but this might not be computationally efficient. Here, we used a faster algorithm (Algorithm 1) to determine the starting and ending time points of admixture.
Let E and S be the ending and starting time points (in generations, prior to the present) of the admixture, which we wanted to search for to minimize the objective function. The search starts from (E 0 , S 0 ) = (1, J), where J is the upper bound for the beginning of the admixture event, which can be set to be a large integer to seek for a relatively ancient admixture event. In our analysis of recent admixed populations, we set J = 500.
where θ 0 , θ 1 can be determined by ordinary least squares. We chose the proposal that resulted in a smaller value for f. The search stopped when the value of f with (E k−1 , S k−1 ) was no larger than that of either proposal or E k = S k . In this way, we could readily estimate the time interval of the admixture event (G end , G start ) quickly.

Algorithm 1:
Result evaluation. To check our assumption of the true history and evaluate the inference, an intuitive way is to compare empirical weighted LD with the fitted LD. Here, we used two quantities: msE and Quasi F, defined by the following: with e i being the ith entry of e. This reflects goodness of fit and strength of background noise. A smaller msE indicates less background noise and better fit.
Scientific RepoRts | 7:43054 | DOI: 10.1038/srep43054 where Ẑ is the fitted weighted LD obtained from iMAAPs, which theoretically can be regarded as the de-noised weighted LD. e ′ is a vector of length I, with the ith entry denoted by ′ e i . We looked at the quasi-F statistic = A reliable result should have both small msE and small F values. Particularly, F is involved in model comparison: when F is too large, one would suspect that the true admixture history is far from any one of these models. Both F and msE are involved in revealing data quality. If F is small but msE is large, one would suspect that the quality of data is not good enough to draw convincing conclusions. Further explanation of these statistics is in Results and Discussion sessions.
Identification of the best-fit model. For the convenience of illustration, we defined the core model as the model used to infer admixture time. When inferring admixture of a target population, HI, GA, CGF1, CGF2, GA-I, CGF1-I and CGF2-I are used as the core models for conducting inference. Because GA-I, CGF1-I and CGF2-I describe more general admixture models than GA, CGF1, and CGF2, we classified model selection into two cases: one case is to identify the best-fit model(s) among the HI, GA, CGF1, and CGF2 models, whereas the more general case is to determine the best-fit model(s) among HI, GA-I, CGF1-I and CGF2-I models. In both cases, the same strategy is adopted, which depends on the pairwise paired difference of pseudo log(msE) values associated with each core model, which will be defined later. For an admixed population, there are N + 1 observed weighted LD curves obtained as follows: N (typically 22) autosomal chromosomes are considered in an individual genome, and one weighted LD curve is calculated from all these N chromosomes while the other N weighted LD curves are obtained by jackknife resampling, leaving out one chromosome for each LD curve 1,10,11 . Next, we fit each observed weighted LD curve for each core model by estimating θ 0 , θ 1 and the time interval, which in turn allowed us to obtain the msE value associated with the optimal parameters for each weighted LD curve. Taken together, a total of N + 1 msE values associated with N + 1 LD curves were evaluated in each core model. For model M, the log(msE) obtained from all N chromosomes was denoted by M 0  and that from the LD curve with the q-th chromosome was left out by   and treated these pseudo values approximately as independent. Next, we defined the best-fit core model(s) to be the model(s) with significantly small ˆq M  . A pairwise Wilcoxon signed-rank test was conducted for the pseudo log(msE) of the four models. More precisely, Wilcoxon signed-rank test was applied to all pairs of models with the ˆq M  being paired by index q, and then the p-values were adjusted to control family-wise error rate (see Table 1). We used the Holm-Bonferroni method to adjust p-values 18 . When  q HI was not significantly larger than those of the best model, i.e., the model associated with the smallest sample median of pseudo log(msE) values, HI was selected because HI is a simpler model compared with the others. Otherwise, the models whose  q M was not significantly larger than those of the best model were selected (the best model was selected as well). The significance level was set to be 0.05. Here, we paired the pseudo values according to index q and used Wilcoxon signed-rank test on the paired differences. ˆq M  is strongly correlated with q and hence q is a major covariate that must be controlled in the test to gain higher power. This is also the reason that even though theoretically there are examples where the best model, according to our definition, can be significantly worse than other models in our process, we still use this method considering that such extreme cases are unlikely in practice. In addition, log(msE) rather than msE was used because after logarithm transformation, the small values of msE could also have huge effect to the comparison. That is to say, we could better detect the   difference between small msE, thus gaining greater power in the test. This claim is also justified by our experience.
In Table 1, we listed the adjusted p-values to determine the best-fit model(s) under various scenarios. In the simulation of HI (100), HI model was inferred as the best-fit model because  q HI is not significantly larger than  - In the cases of CGF1 (1-50), GA-I, CGF1-I, and CGF2-I were inferred as the set of best-fit models because we cannot distinguish the best fit model from GA-I, CGF1-I, and CGF2-I models. This case was marked as "Undetermined" in Table 2.

Software
Our algorithm has been implemented in an R package 19 , named CAMer (Continuous Admixture Modeler). The package is available on the website of population genetic group: http://www.picb.ac.cn/PGG/resource.php or on Github: https://www.github.com/david940408/CAMer.

Results
Simulation studies. Admixed populations were simulated in a forward-time way under different admixture models with the software AdmixSim 20 , which is under the framework of copying model that new haplotypes are assembled from the segments of the source populations' haplotypes generation by generation 4,21 , and the same simulation strategy has been used in the previous work 4 . Simulation was initiated with the haplotypes from source populations (YRI and CEU) and the haplotypes for the admixed population were generated by resampling haplotypes with recombination from source populations and the admixed population of last generation. During the simulation, population size was kept as 5000 and migration rate was controlled by the admixture model with the final admixture proportion in the admixed population to be 0.3. We employed a uniform recombination map in our simulation, which means recombination rate between two markers is positively proportional to their physical distance. For each model, simulation was performed using 10 replicates; each replicate contained 10 chromosomes with a total length of 3 Morgans. To evaluate the performance of our algorithm, we simulated admixed populations under the following conditions: (1) HI of 50 and 100 generations, designated as HI (50) and HI (100), (2) GA of 50 and 100 generations, designated as GA (1-50) and GA (1-100), respectively, (3) CGF of 50 and 100 generation, population 1 as the recipient, designated as CGF1 (1-50) and CGF1 (1-100) respectively, (4) CGF-I of a 70-generation admixture followed by 30-generation isolation, and a 30-generation admixture followed by a 70-generation isolation, with population 1 as the recipient, designated as CGF1-I (30-100) and CGF1-I (70-100) respectively, and, (5) GA-I of a 70-generation admixture followed by a 30-generation isolation and a 30-generation admixture followed by a 70-generation isolation, designated as GA-I (30-100) and GA-I (70-100), respectively.
With simulated admixed populations, we first used the HI, GA and CGF models as core models to conduct inference (see Supplementary Fig. S1). When the simulated model was a HI, GA, or CGF model, our method was able to accurately estimate the admixture time, as well as to determine the correct model, with an accuracy of 73.33%. When the simulated model was a CGF-I or GA-I model, the estimated time based on the core model HI was within the time interval of the admixture, whereas all best-fit models were HI (see Table 2).
With the same set of simulated admixed populations, we also used AdmixInfer 9 to determine the admixture model and estimate admixture time, which is based on the length distribution of CAT. To avoid any errors introduced by haplotype phasing and local ancestry inference, we analyzed the ancestral segments generated from AdmixSim. We found that AdmixInfer attained pretty accuracy in determining the admixture model and estimating admixture time when the simulation is under HI, GA, or CGF model. However, it could only give HI model as the best-fit model when the simulated admixture is under GA-I or CGF-I model. (see Supplementary Table S2) These results indicated the limitation of using the GA and CGF models in inferring admixture history, no matter the information from LD or CAT is used for inference.
We next employed GA-I, CGF-I and HI as core models for performing inference (see Fig. 3 and Supplementary Figs S2-11). With HI, GA, or CGF being considered as the true model, our estimation of the optimal model remained accurate. On the other hand, when the true model was GA-I or CGF-I, the failure rate decreased by 35%, compared to the estimation in the previous setting, but it was still at a very high level.  Table 2. Accuracy of model determination. Here, as our method can hardly distinguish CGF1 from CGF2 model, we regard CGF1, CGF2 as the CGF model; CGF1-I and CGF2-I as the CGF-I model, which are different from GA-I and HI models. Here, "correct" denotes the best-fit model is the true model; "Undetermined" means the true model can not be determined from the best-fit models; "Wrong" denotes the true model is not given. Furthermore, the estimated time intervals were wider than those of the true ones, although the results were still more accurate than those using GA and CGF as core models (see Table 2). By introducing the GA-I and CGF-I models as core models, CAMer can resolve the admixture into continuous time interval. Considering that CAMer is not so powerful in determining the best-fit admixture model (see Tables 1 and 2), in empirical studies, we presented the results from CAMer with estimations on all core models and the model(s) fitting best the data.

Empirical analysis.
We applied CAMer to the selected admixed populations from HapMap, HGDP, and 1KG. For each target population, we first used iMAAPs to calculate the weighted LD and fit the weighted LD decay curve with a numeric method 11 . Next, with the weighted LD of target populations, we determined the admixture model and estimated admixture time with CAMer. Quasi F and msE are designed for evaluating the inference with CAMer. The value of msE usually indicates data quality: small msE may indicate a high signal-to-noise ratio (SNR) and vice versa. The quasi F value measures the goodness of fit of the model we employed to fit the admixture event. A small F value indicates that the model we used was of satisfactory performance in modeling an admixture event. In our analysis, we used 10 −5 as the threshold for msE and 1.5 for F. Therefore, when the msE value ≤ 10 −5 and the F value ≤ 1.5, we could not "reject the null hypothesis" that the related model was the true model, i.e., the model well fit the admixture event. On the other hand, an msE value ≥ 10 −5 indicates low-quality data that is incapable of identifying the best-fit model, whereas a F value ≥ 1.5 prompts us to "reject the null hypothesis" and concludes that the model does not well fit the admixture. In the case of the same population from different databases, the data with smaller msE values were given more credits. For example, we obtained samples of ASW from the HapMap and the 1KG. With the ASW data (CEU and YRI as source populations) from HapMap, the best-fit model was GA-I of 2-8 generations, and both msE and F values indicated that the inference was acceptable (see Supplementary Fig. S12). Similarly, using the ASW data (CEU and YRI as source populations) from 1KG, the best-fit model failed to be determined among GA-I, CGF1-I, and CGF2-I (see Supplementary Fig. S13). However, all the quasi F values bigger than 1.5 indicated that these models did not satisfactorily fit the admixture event. Because the msE value of the data set from 1KG was smaller, the conclusion using ASW was as follows: based on the best data we had, the time intervals estimated under the HI, GA-I, CGF1-I, and CGF2-I model were 5 generations, 2-9 generations, 1-11 generations, and 1-9 generations, respectively. Furthermore, none of these models satisfactorily modeled the admixture, whereas the HI model showed better performance. We also applied CAMer to other admixed populations (see Table 3

Discussion
Modeling the demographic history of an admixed population and estimating time points of admixture event are essential components of evolutionary and medical research studies [5][6][7][8][9]11,22 . Previous methods have employed the length distribution of ancestral tracts [6][7][8] , which highly depends on the accuracy of local ancestral inference and haplotype phasing. Another limitation is that only HI, GA, and CGF models were utilized to fit the admixture as well as in identifying the best-fit model. In the present study, our simulations showed that when the true model was not HI, GA, or CGF, the generated inferences were relatively difficult to interpret.
Our method, CAMer, can be utilized in inferring admixture histories based on weighted LD, which can be calculated using genotype data with iMAAPs 11 . Furthermore, we extended the GA and CGF models to the GA-I and CGF-I models in order to infer the time interval for a period of continuous admixture events followed by  isolation. Although HI model is a degenerate case for both GA-I and CGF-I models, where the admixture window becomes 1 generation, we kept it in our method because it is the most popular model employed in previous admixture studies. Considering the difficulty in fitting problem with polynomial functions, it is in our expectation that CAMer was not consistently accurate in determining the admixture model based on the weighted LD decay. However, its natural advantage of independence of both haplotype phasing and local ancestry inference makes it privilege to other CAT based methods. Our simulations indicated that its time interval estimations were reliable when its assumption that the true admixture history could be well approximated by one of the core models is valid. Two quantities, namely msE and quasi F, were used to check the assumption of our method stated above and evaluate the credibility of the models' inference. These two quantities should both be taken into consideration to determine whether the models well described the admixture history. Both the data quality and the goodness of fitting of models can affect the value of msE, although the F value mainly measures the goodness of modeling. Informally, for the convenience of interpretation, msE can be an indicator of data quality, while F value can be used to check model assumption on admixture history. In our analysis, we suggested thresholds for msE and F to determine whether the null hypothesis should be rejected or not, which may be too strict in empirical analysis. Actually, msE and F values together measure whether the observed weighted LD can be well fit by the best-fit model(s). For example, the fitting process showed poor performance in the MKK population, which was accompanied by exaggerated msE and F values, showing significant inconsistencies between the observed and fitted weight LD curves, which indicates that the true admixture history cannot be well explained by any of the core models (see Supplementary Fig. S17). Therefore, in empirical analysis, it can be informally considered that the msE value reflects the quality of the data, whereas F value describes the performance of the model, although both of them measure the goodness of fitting.
In our previous study 11 , we fit the weighted LD with high degree polynomial functions. However, this approach did not fully reveal the occurrence of continuous admixture. To address this issue, the present study developed CAMer to model admixture as a continuous process. CAMer also employed extensions of the classic continuous models, GA-I and CGF-I, which may bring the bias to have a wider admixture window when the real admixture exists in a short time. As we discussed earlier, another limitation for CAMer is its poor performance to determine the correct admixture model. Therefore, in empirical data analysis, we suggest all core models, rather than the best-fit model(s), should be examined. Taken together, despite there is space to further improve in the future, CAMer is a powerful method to model a continuous population admixture, which in turn would help us elucidate the complex demographic history of population admixture.