## Abstract

Migrations have played an important role in shaping the genetic diversity of human populations. Understanding genomic data thus requires careful modeling of historical gene flow. Here we consider the effect of relatively recent population structure and gene flow and interpret genomes of individuals that have ancestry from multiple source populations as mosaics of segments originating from each population. This article describes general and tractable models for local ancestry patterns with a focus on the length distribution of continuous ancestry tracts and the variance in total ancestry proportions among individuals. The models offer improved agreement with Wright–Fisher simulation data when compared to the state-of-the art and can be used to infer time-dependent migration rates from multiple populations. Considering HapMap African-American (ASW) data, we find that a model with two distinct phases of “European” gene flow significantly improves the modeling of both tract lengths and ancestry variances.

DNA sequencing is an invaluable tool for understanding demographic relationships between populations. Even with a limited number of genetic markers, measured across individuals and populations, it is often possible to estimate relatedness between populations, ancestry proportions in admixed populations, or sex-biased gene flow. The availability of dense genotyping platforms and high-throughput sequencing technology has enabled refined analyses of genetic diversity.

Because of recombination, different loci along an individual genome can reveal different aspects of its ancestry. Consider a sample and its ancestral population at some time *T* in the past, and suppose that we give ancestral individuals subpopulation labels, defining *source* populations. These labels are typically chosen to represent subgroups that have increased genetic homogeneity due to cultural or geographic reasons. Then a simple summary of the demographic trajectory of a sampled allele is the source population from which it originated. We say that an individual is “admixed” if it draws ancestry from multiple source populations—thus admixture is not an intrinsic property of individuals, but depends on our choice of labels and time *T*. An example of subpopulation labels often used to study human populations in the Americas are the European, Native American, and West African populations prior to the advent of massive intercontinental travel. Many routines have been proposed to infer the source population along the genome of admixed individuals (Ungerer *et al.* 1998; Tang *et al.* 2006; Falush *et al.* 2003; Hoggart *et al.* 2004; Patterson *et al.* 2004; Sankararaman *et al.* 2008; Bercovici and Geiger 2009; Price *et al.* 2009). These typically proceed by locally matching an admixed genome to panel populations chosen as proxies for the source populations, revealing a mosaic of tracts of continuous ancestry (Figure 1). In this work we use PCAdmix (Brisbin 2010), a heuristic approach for local ancestry inference. PCAdmix first divides the genome in windows of typical width of 10–50 kb. For each window, the probability that the sample haplotype originates from any of the panel populations is estimated on the basis of the position in PCA space. Finally, PCAdmix uses these probabilities as emission probabilities of a hidden Markov model and ancestry is inferred via Viterbi decoding.

Local ancestry patterns have been used to identify disease loci (see Seldin *et al.* 2011 and references therein) and to search for regions experiencing selection (Tang *et al.* 2007; Bhatia *et al.* 2011). They also provide hints about the history of migration (Pool and Nielsen 2009). The purpose of this article is to understand and model the observed ancestry patterns on the basis of detailed demographic models, to learn about human demography, and to empower selection and association scans. In particular, we are interested in the length distribution of the continuous ancestry tracts, and the variation in ancestry proportions across chromosomes and individuals.

A dominant stochastic process leading to these patterns is recombination, which tends to break down segments of continuous ancestry in admixed individuals. As a result, the length of continuous ancestry tracts tends to be shorter for more ancient admixture. The tract-length distribution is sensitive to details of recent migration (*i.e.*, tens of generations) and is thus complementary to analysis based on the joint site-frequency spectrum (Gutenkunst *et al.* 2009; Gravel *et al.* 2011), which is more sensitive at time scales of hundreds to thousands of generations.

Recently, Pool and Nielsen (2009) proposed a model in which a target population receives migrants from a source population, initially at a constant rate *m*_{2}. Starting at a time *T* in the past, the rate changes to *m*_{1}. In this model, back migrations are not allowed, recombinations within migrant chromosomes are neglected, and tracts shorter than a cutoff value are forgotten (since migration occurs over an infinite period, this is necessary to avoid having a genome completely replaced by migrants). Assuming that recombinations occur according to a Poisson process, these approximations allow for an analytical solution for the distribution of tract lengths, which was used to infer demographic events in mice (Pool and Nielsen 2009). This model is limited to admixture proportions weak enough that recombinations between migrant chromosomes are unlikely. A second limitation is that the model assumes two epochs of constant migration rate, which might or might not be the most appropriate for a given population. The special case *m*_{2} = 0 has been used to infer demographic histories in humans for North African individuals (Henn *et al.* 2012).

Here we propose a more general approach to predicting the distribution of tract lengths that can accommodate both time-dependent and strong migration. This approach builds on that of Pool and Nielsen (2009) but introduces multiple improvements. First, general time-dependent migrations can be considered. Second, recombinations between tracts of the same ancestry are not neglected, allowing for the modeling of strong migration and the simultaneous study of tracts of multiple ancestries. Third, chromosomal end effects are modeled explicitly. Fourth, our model can be modified to incorporate errors in tract assignments. As in the Pool and Nielsen approach, we model recombination as a Poisson process with a unit rate per morgan, and the recombination map is taken to be identical across populations [a reasonable approximation at the centimorgan scale (Wegmann *et al.* 2011)]. To perform demographic inference, we further require that local ancestry inference can be performed to high accuracy using one of the methods mentioned above. Whether this can be done depends on the degree of divergence of the ancestral populations (or sources), the availability of data for panel populations that are good proxies for the sources, and the possibility of accurately phasing diploid genomes.

Admixture history also leaves a trace in the variance in admixture proportions across individuals, as stochastic mating and recombination tend to uniformize ancestry proportions with time (Verdu and Rosenberg 2011). Generalizing the models of (Ewens and Spielman 1995; Verdu and Rosenberg 2011) to include the effects of recombination in a finite genome and drift, we show that after a discrete admixture event, the variance decays in time in three consecutive regimes, first exponentially as differences in individual genealogies average out, then linearly as recombination creates shorter tracts, and finally exponentially again as drift fixes local ancestry along a chromosome. A simple approximate equation captures all three regimes accurately. By contrast, variance in continuous migration models is dominated by the first regime, and the expressions from the model of Verdu and Rosenberg (2011) are reasonably accurate (see Appendix 3).

In general, distinguishing the effects of population structure and time-dependent patterns of gene flow is not straightforward, and the inference problem is prone to overfitting, as is the case, *e.g.*, for inference based on the site-frequency spectrum (Myers *et al.* 2008). However, our analysis shows that tract lengths, and more generally ancestry correlation patterns, can help resolve subtle differences in patterns of historical gene flow. An implementation of the proposed methods for tract length modeling, called *tracts*, is available at http://tracts.googlecode.com.

## Theory

### Admixture models: definitions and global properties

We construct a model for the admixture of diploid individuals that takes into account recombination, drift, migration, and finite chromosome length. Since a full coalescent treatment of these effects is computationally prohibitive (Griffiths and Marjoram 1996), we simplify the model to consider only the demography of our samples up to the first migration event, *T* generations ago. We label generations *s* ∈ {0, 1, 2, …, *T* − 1}, and the total fraction of the population *m*(*s*) that is replaced by migrants in a generation *s* can be subdivided in contributions *m _{p}*(

*s*) from

*M*migrant populations:

*p*∈ {1, …,

*M*}. We treat the replacement fraction

*m*(

_{p}*t*) as deterministic, while the replaced individuals are selected at random (see Figure 2). Generations follow a Wright–Fisher model with random mating in a population with 2

*N*genome copies, each with

*K*finite chromosomes of morgan length {

*L*}

_{i}

_{i}_{=1…,}

*. We consider two different variations of the Wright–Fisher model with recombination.*

_{K}The first variation (model 1) is meant to be the most biologically motivated and will be used for all simulations. Starting from a finite parental diploid population of size *N*, we first replace *m*(*s*)*N* randomly selected individuals with diploid migrants. Diploid offspring are generated by drawing one gamete from each of two randomly selected diploid parents. Gamete formation is a Markov path with transition rate of one transition per morgan across the two parental chromosomes (see Figure 2B).

Model 1 results in long-range, non-Markovian correlations along the genome. This complicates the modeling without necessarily having a large effect on most global statistics. We therefore also consider a more tractable model (model 2) in which gametes are drawn from the migrant populations with probability *m*(*s*) and are otherwise generated by following a Markov path along all nonmigrant parental gametes (see Figure 2). The reason for singling out new migrants is that it is possible to generate their gametes as in the more realistic model 1, without sacrificing tractability. Model 2 may not capture all long-range correlations in ancestry but it has the correct distribution of crossovers and for small portions of the chromosomes is very similar to model 1: the only difference is that each draw from the parental gamete pool is independent in model 2, whereas the fact that a diploid individual can have multiple offspring induces a small degree of correlation between draws in model 1. Unless otherwise stated, we calculate all population-wide statistics after the migration step, but before gamete generation.

Model 2 is reminiscent of the Li and Stephens (2003) copying model used in HAPMIX (Price *et al.* 2009), as it also neglects back-and-forth recombinations due to multiple crossovers during a single meiosis. The purpose of the models are different, in that the current Markov models attempt to simulate gamete formation from parental chromosomes and represent evolution in time, whereas the Li and Stephens model attempts to simulate an unobserved haplotype on the basis of haplotypes from the same generation. The Markov ancestry transition model used in HAPMIX (and many other local ancestry inference software) corresponds to a special case of model 2 when each population contributes migrants at a single generation.

Local ancestry patterns are sensitive to the three stochastic processes of migration, recombination, and random genetic drift. Where possible, we take all three effects into account. By contrast, we do not model the effects of population structure, of selection, and of population size fluctuations. We derive our results under the assumption that local ancestries can be determined exactly; the effects of misidentification are discussed throughout, together with possible correction strategies.

Given a history of migrations, it is relatively straightforward to calculate the expected population averages for ancestry proportions and tract lengths. If *m*(*s*) is the total fraction of the population that is replaced by migrants, *s* generations ago, with *m _{i}*(

*s*) from population

*i*, the expected ancestry from population

*i*at a time

*t*in the past is the sum over generations

*s*of migrant contributions

*m*(

_{i}*s*) weighted by the survival probability

*t*. After the migration step, the ancestry proportions are

We can follow a similar procedure to obtain the expected density *w _{ij}* of ancestry switch points from population

*i*to population

*j*per morgan, replacing the amount of new migrants

*m*(

_{i}*s*) by the density of new switch points, which are proportional to the recombination rate (assumed constant with unit rate in genetic units) and the expected fraction of the genome

*h*(

_{ij}*s*) that is heterozygous with respect to ancestries

*i*and

*j*after generation

*s*. In the gamete pool, we findThe ancestry heterozygocity

*h*can be evaluated using a recursive equation (such as Equation A1), as in the case of allelic heterozygocity. In the absence of drift,

_{ij}*h*(

_{ij}*s*) = (1 −

*m*(

*s*))

*α*(

_{i}*s*+ 1)

*α*(

_{j}*s*+ 1). In the population (before gamete generation), the sum over

*s*starts at

*t*+ 1 rather than

*t*. The expected number of switch points per morgan at time 0 is therefore

To estimate the expected tract length for ancestry *i* on a chromosome of length *L*, we divide the expected length covered by this ancestry, *α _{i}*(0)*

*L*, by the expected number of tracts of this ancestry, which is

### Tract-length distribution

For illustration, we first consider a source population (Blue), and a target population (Red), with a single, infinitely long diploid chromosome. At generation *t* = *T* − 1, a fraction *m* of population Red is replaced by individuals from population Blue. Consider the Markovian Wright–Fisher model discussed above (model 2). In this model, the position of the closest recombination to either side of a point along an infinite chromosome is exponentially distributed and there is no memory of previously visited states along a chromosome. The chromosomes resulting from this admixture process can therefore be modeled as a continuous-time Markov model with a Red and a Blue state (Figure 3A), where each recombination event corresponds to a Markov transition and the continuous Markov “time” corresponds to the position along the chromosome. The transition rate out of a state in this model is proportional to the number of recombinations, namely *t* − 1 per morgan: since recombinations within first-generation migrants do not induce ancestry changes, and we suppose that we sequence somatic cells at generation 0, recombination can occur only during gamete formation at generations 1, ..., *t* − 1. If a recombination occurs, the probability of transitioning is *m* to the Red state and is (1 − *m*) to the Blue state.

We are interested in the length distribution of continuous segments in the Blue or Red ancestry, independent of the number of within-ancestry transitions, which are difficult to detect. We avoid these complications by setting the self-transition rates to zero: this does not affect the trajectories, but now all transitions change the ancestry. We therefore have the model shown in Figure 3A, and the distribution of tract lengths *φ _{i}*(

*x*) is equal to the exponential distribution of distance between Markov transitions:

*t*= 1, since this situation produces tracts that are infinite in the infinite-chromosome limit.

#### Multiple populations, discrete migration:

As long as the migration from each population is limited to a single generation and the target population is infinitely large, model 2 produces Markovian trajectories along ancestry states. To see this, consider a point *x* along the genome in a segment from ancestry *p* that arrived *t* generations ago. As before, the distance to the first recombination event downstream from *x* is exponentially distributed (with rate *t* − 1), and the timing *τ* of the recombination is uniform on (1, *t* − 1). Moreover, since gametes in model 2 are formed by following a Markov path in the parental gamete pool, the probability of observing ancestry *p*′ downstream from the recombination is proportional to the ancestry proportions in the parental pool at the time *τ* of the recombination. Thus we have the discrete transition rate*t* of ancestry *p*. We note that the Markov property over ancestry states would be lost in model 1, because the state downstream of the recombination is correlated with upstream states. Drift reduces the transition rates and also breaks the Markov property: mitigation strategies are discussed in Appendix 1. The Markov property over ancestry states is also lost if a population contributes migrants over many generations, and our next step is to restore the Markov property in this situation by extending the state space.

#### General incoming migration in the absence of drift:

We now allow for general incoming migration histories that start at a time *T* − 1 in the past. For each generation *t* ∈ {0, …, *T* − 1}, a fraction *m _{p}*(

*t*) of the individuals from the target populations are replaced by individuals from the source population

*p*, with

*m*(

*T*− 1) = 1. Since the ancestry switches are no longer Markovian in the general migration case, it is convenient to consider states defined by both ancestry

*p*and time of arrival

*t*. Intuitively, we may imagine that we have a large number of migrant populations (

*p*,

*t*), each contributing migrants over a single generation (see Figure 3, B and C). Here the Markov property is maintained, but ancestry states can now correspond to multiple Markov states.

We first calculate the transition rates between states (*p*, *t*) as we did for the discrete migration case. First, the probability of encountering state (*p*, *t*) downstream from a recombination that occurred at time *τ* is

As before, given a point *x* in state (*p*, *t*), the position of the next downstream recombination is exponentially distributed with rate *t* − 1, and the time of this recombination is uniformally distributed on (1, *t* − 1). In the two Wright–Fisher models considered here, states on either side of the recombination are uncorrelated, and we can write the discrete transition probabilities*p*. The continuous transition rate is obtained by multiplying the discrete transition rate by the continuous overall transition rate *t* − 1:

Given the transition matrix *Q*, we can use standard tools for the study of Markov chains to efficiently estimate the length distribution of excursions on Markov states corresponding to a single ancestry. In Appendix 2, we first derive results under the approximation that chromosomes are infinitely long. We account for finite chromosomes by studying the distribution of tract lengths in finite windows, randomly chosen along the infinite chromosomes. We thus obtain a distribution of tracts *φ _{p}*(

*x*) for each population

*p*. To compare these predictions to observed data, a computationally efficient strategy is to bin data by tract length and treat the observed counts in each bin as an independent Poisson variable with mean obtained by integrating

*φ*(

_{p}*x*) over the bin range.

Short ancestry tracts are likely to have both elevated false-positive and false-negative rates, and inference based on such tracts is likely to be biased, whereas longer tracts can be detected with increased confidence. Following Pool and Nielsen (2009), we therefore perform inferences using only tracts longer than a cutoff value *C*. We should emphasize that a large number of uniformly distributed spurious short tracts may still affect the distribution of longer tracts, making nonexponential distributions look more exponential. Therefore, significant assignment error may cause an underestimation of the amount of continuous migration. By contrast, drift would tend to reduce the transition rates and cause underestimates of the time since admixture (see Appendix 1).

### Variance among individuals

We now consider the variance across individuals in total migrant ancestry *X ^{p}* from population

*p*, measured as a proportion of the morgan length of the genome whose origin is from

*p*. The variance in ancestry can be separated in two components, which we label the genealogy variance and assortment variance. The genealogy variance is due to a different number of migrant ancestors; if a randomly chosen fraction

*m*of the population is replaced by migrants at each generation, a fraction

*m*

^{2}of individuals will have two migrant parents, 2

*m*(1 −

*m*) will have one migrant parent, and (1 −

*m*)

^{2}will have none. The assortment variance accounts for the fact that two individuals with the same genealogy can vary in their genetic ancestry proportions, since not all ancestors contribute the same amount of genetic material to an individual. Recombination and the independent assortment of chromosomes tend to reduce such variance.

We can use the law of total variance, conditioning over the genealogies *g*, to isolate these two contributions to the variance Var(*X ^{p}*):Here is the fraction of migrant ancestry from population

*p*, based on the genealogy

*g*. Alternatively, this can be thought of as the infinite-sites expectation for the ancestry proportions. The first term therefore represents the genealogy variance in ancestry, whereas the second term represents the assortment variance. Because of random chromosome assortment, the variance in ancestry among chromosomes is informative of the assortment variance. We discuss in Appendix 3 how, in the absence of drift, the variance can thus be broken down in these two components without requiring a demographic model. We discuss below how to obtain expectations for each components given a specific demographic model.

#### Genealogy variance:

To ease calculations of the genealogy variance, we neglect correlations due to overlap between individual genealogies and describe each individual as being sampled from an independent genealogy (in a randomly mating population, this amounts to neglecting drift). In this model, the genealogy variance Var* _{g}* is easily calculated. Considering the genealogy

*g*of a nonmigrant sample up to

*T*generations ago (we label the current generation 0, and the generation with the first migrants

*T*− 1), we first note thatwhere

*i*and 0 otherwise. Results with continuous admixture since time immemorial can be obtained by taking a limit

*T*→ ∞. In such cases, the approximation of independent pedigrees eventually breaks down, but the resulting expression might remain approximately correct if the majority of present-day genomes originate from recent migrants.

The expectation over genealogies *g* and assortments is then *α _{p}*(0). The calculation of is also straightforward if we can calculate the expectation . For

*i*and

*j*, or one migrant on each of the separate branches,with

*d*half the tree distance between leafs

_{ij}*i*and

*j*. Then we can write the sum over half distances, weighted by the number of leaf pairs at each distance:Since , we haveIn the two-population pulse model, with

*m*

_{p}_{=1}(

*t*) =

*mδ*

_{t}_{,}

_{T}_{−1}, we have the expected with a rapid exponential decay of the variance as a function of

*T*. By contrast, if we have continuous migration of population

*p*in a target population, with,

*T*. Finally, in the case in which two populations provide respectively

*pm*and (1 −

*p*)

*m*migrants to a target population at each generation since the beginning of time, we have the simple expressionThis expression supposes that the variance is calculated after migration occurs. If variance is calculated before replacement by migrants, the factor of 2 disappears, and we recover Equation 47 in Verdu and Rosenberg (2011).

#### Assortment variance:

To study the global ancestry variance due to assortment, a natural starting point is to consider the ancestry variance at a particular point in the genome. In a randomly mating population with two ancestries, the variance in ancestry at a site is *h*/2, where *h* is the ancestry heterozygocity at that site. The ancestry heterozygocity can be calculated using the same recursive strategy commonly used for allelic heterozygocity (Equation A1). The case of three or more ancestries can be reduced to two ancestries by singling out one ancestry and pooling the others. As a specific example, in the case of a pulse migration with migration rates *m* and 1 − *m* at generation *T* − 1, the heterozygocity at generation 0 is*n* tracts by uniformly drawing *n* − 1 separators. We suppose that the ancestry is chosen independently on each segment, with variance *h*_{0}/2. Then the variance in ancestry in the large-*n* limit isThe effect of drift is therefore captured by the decay of ancestry heterozygocity with time, whereas the effect of recombination is captured by the number of independent tracts *n*, which is proportional to the number of recombinations. In the case of a pulse of migration *T* generations ago without drift, we write *n* = 1 + (*T* − 2)*L _{i}* for a single haploid chromosome (the 1 accounts for the chromosome edge and can be neglected for large

*TL*), and 2

_{i}*K*+ 2(

*T*− 2)

*L*for a diploid genome with

*K*chromosome pairs of total length

*n*, this expression provides excellent quantitative agreement with simulations over multiple regimes (Figure 6). If we model the variation over time of the population ancestry proportion as a random walk with decreasing step size

*σ*= 0.5 ± 0.1, a relatively modest uncertainty given the small population size. Assortment variance for continuous migration models is discussed in Appendix 3.

## Comparison with Simulation and Experimental Data

In this section we first present results of Wright–Fisher simulations, comparing our model predictions to the simulation results. We then consider the HapMap African-American panel, for which we performed local ancestry inference and analyzed the tract length distribution.

### Tract lengths

We performed a 30-generation diploid Wright–Fisher simulation (using model 1; see Figure 2) of 10000 chromosomes of length 1 M with continuous gene flow from population 1 into a population initially composed of individuals from population 2. We considered three different migration intensities, namely *m*_{1} = 0.001, 0.03, and 0.05 per generation. We kept track of the ancestry of each segment during the simulation, so that the continuous ancestry tracts could easily be tabulated. In Figure 4, we compare the observed histograms of tract lengths for population 1 (dots) to predictions from Equation 10 in Pool and Nielsen (2009) (dashed lines) and to predictions from the Markovian Wright–Fisher model (model 2 on Figure 2), using rates from Equation 2 and implemented as described in Appendix 2 to account for finite chromosome length (solid lines). As expected, the predictions of the two models are similar when migration rates are low and differ substantially when we depart from the weak migration assumptions of the Pool and Nielsen model (see Figure 4). The Markov model predictions are in good agreement with the simulations over the range of models considered, including when the migrant population becomes the majority population.

We now consider the HapMap African-American panel (ASW) (International HapMap 3 Consortium 2010) and focus on 20 unrelated samples that were trio phased, to reduce biases due to phasing errors. We obtained local ancestry inferences using PCAdmix (Brisbin 2010), using 132 unrelated HapMap samples from Europe (CEU) and 204 from West Africa (YRI) as reference panels. We used windows of size 0.3 cM for the HMM and based our inferences on the number of tracts longer than 10 cM. We pooled tracts in 50 bins according to tract length (chromosomes with no ancestry switches were in a separate bin independent of the chromosome length) and calculated model likelihood assuming that counts in each bin are Poisson distributed with mean given by the model predictions for this bin.

We compared inferences on the basis of two different models; (a) a “pulse” model, with a single migration event, and (b) a two-pulse model, with a subsequent migration of Europeans (Figure 5). Model b has two additional parameters, corresponding to time and proportion of the subsequent European migration. A likelihood-ratio test shows that ln(ℒ_{b}/ℒ_{a}) ≃ 7. To establish the significance of the extra two parameters, we simulated 1000 random tract-length distributions from the maximum-likelihood model a, and obtained maximum-likelihood estimates for both models. The probability of obtaining such a likelihood ratio under model a is *P* = 0.002.

### Ancestry proportions and variance

Simulations of 80 individuals, each with 22 autosomal chromosomes of realistic lengths (namely 2.78, 2.63, 2.24, 2.13, 2.04, 1.93, 1.87, 1.70, 1.68, 1.79, 1.59, 1.73, 1.27, 1.16, 1.26, 1.35, 1.30, 1.19, 1.08, 1.08, 0.62, 0.73 M, for chromosomes 1–22, respectively) and 30% of initial admixture proportion, illustrate many of the effects predicted in variance models. The global ancestry proportions and fraction of sites heterozygous for ancestry fluctuate considerably over the first few generations, but the fluctuations decrease in time as ancestry proportions approach a fixed value and ancestry heterozygocity decays following Equation 7.

Figure 6 shows that the variance in ancestry across individuals follows three different regimes; first, the variance is dominated by the genealogy variance, with a rapid exponential decay. After about 10 generations, the assortment variance starts to dominate, and decays polynomially due to recombination until drift becomes important, where an exponential decay is resumed, although at a much reduced rate.

Equation 8 captures these three regimes in quantitative detail—the average variance over 50 independent simulations follows the model prediction closely. The continuous migration case, in which genealogy variance tends to dominate, is discussed in Appendix 3.

Comparing the ancestry variance from the African-American data to those predicted by the demographic models, we find that the pulse model predicts a genealogy variance of 0.0005, whereas the variance in the model with two distinct pulses is 0.002. The total variance in the African-American sample is 0.0047, of which we infer that 0.0041 is due to genealogy variance (using the method described in Appendix 3). Thus the model with two pulses of migration is again more realistic than the single pulse model; the fact that it still underestimates the variance can be due to a combination of factors that have not been modeled: our demographic model may be underestimating low level, very recent migration because of the parameterization as two discrete pulses of migration, and both population structure and errors in ancestry assignment may be adding to the observed variance.

## Discussion

### Limitations and possible improvements

A limitation of all demographic inference methods is that the model space is very large, and the information available to learn about the models is limited. Thus we need to coarsely parameterize model space at the risk of introducing biases. This is similar to the modeling of allele frequency distributions: even though the vast majority of scenarios are inconsistent with the data, the number of models that are consistent with the data remains large, and model fitting often requires simplifying assumptions. When applied to HapMap trio-phased African-American data (ASW), inferred parameters were reasonable and we found evidence for migration patterns that depart from the migration pulse with subsequent random mating, which is at the heart of many approaches. However, distinguishing between continuous migration and nonrandom mating remains challenging.

The demographic inference strategy we presented requires accurate local ancestry assignments. Since longer tracts contain more ancestry information, we expect the most significant types of misassignment to be short, spurious ancestry tracts and the failure to identify real, short ancestry tracts. In the HapMap ASW data used here, the source populations are diverged enough that assignments are relatively reliable down to relatively short tracts, and indeed we find that the number of tracts predicted by the model is in good agreement with the data for the shortest tracts, even though these were not used in the fitting procedure. If the expected number of misidentified short tracts is large enough that it will strongly affect the distribution of longer tracts (by introducing spurious breaks in longer tracts), the Markov models should be modified to include misidentification states, and transition rates could be estimated via simulated admixed individuals.

Alternatively, it is possible to circumvent the local ancestry inference step altogether by focusing on a derived statistic, such as the decay of correlation in ancestry informative markers with genetic distance. Such a method was proposed in Reich *et al.* (2009), for the case of pairwise ancestry correlations in a pulse migration model. Even though such an approach avoids possible biases due to local ancestry assignment, pairwise ancestry correlations become noisy as distance is increased and are thus less sensitive to continuous gene flow patterns. The Markov models presented here provide a natural framework to generalize linkage-based models for more general admixture scenarios, as arbitrary order linkage statistics can be derived in the Markov framework. Furthermore, HMM approaches could be developed to model the complete mosaic pattern without the need to focus on summary statistics such as the tract-length distributions and ancestry variances. Even though such approaches would be more computationally intensive, they may increase the accuracy of the inference, especially when assignment errors are important.

### Conclusion

Overall, we found that the proposed models accurately describe the distribution of ancestry tract lengths and variances when compared to Wright–Fisher simulations. The models we used allow for general migration histories, yet are tractable and can be used for inferring demographic parameters in real data. They are therefore useful to improve our understanding of the consequences of gene flow and our ability to infer demography in populations with complex histories. Such populations have often been underrepresented in medical genetic studies, in part because of complications in the modeling of genetic heterogeneity. As medical genetics sampling efforts strive to reduce this disparity, detailed models for genetic diversity will be increasingly important to make the most out of the resulting data.

## Acknowledgments

I thank Carlos D. Bustamante, Jake K. Byrnes, Brenna M. Henn, Jeffrey M. Kidd, and Damien Simon for useful discussions. This publication was made possible by National Institutes of Health (NIH)/National institute of General Medical Sciences grant number 1 R01 GM090087-01 and NIH/National Human Genome Research Institute grant number U01HG005715.

## Appendix 1: The Effect of Drift on Ancestry Transitions

Drift increases the probability that recombinations occur between segments of the same ancestry. In the infinite-time limit, ancestry will be fixed at every site, no more ancestry switches are created, and the tract-length distribution is constant in time. In the presence of drift, the ancestry switches are no longer Markovian; if a recombination occurs between two IBD segments, it increases the posterior probability that the next recombination will also be between IBD segments. However, it is likely that a Markovian approximation will remain accurate for moderate drift if we take into account the reduced probability of ancestry-switching recombinations.

We first wish to obtain the fraction of recombinations that occur within segments (*p*, *t*), of ancestry *p* having migrated at generation *t*, as these recombinations do not induce ancestry switches and will be most affected by drift. In other words, we want to find the fraction of sites that are homozygous for the ancestry (*p*, *t*), and contrast this to the case with no drift. For this purpose, we consider all other ancestries as a single allele, and in the first step we compute the total homozygocity of nonmigrants in this system *s* generations ago:

This recursion can be initiated with the homozygocity one generation after *t*, namely *c _{p}*

_{,}

*of nonmigrant sites that are homozygous for the*

_{t}*p*,

*t*ancestry at generation

*s*, we write

*c*

_{p}_{,}

_{t}In the driftless case, the probability of the state to the right of a recombination depended only on the time of the recombination. Due to the possibility of recombining within segments identical-by-descent, this is no longer the case when drift is present. However, consider a given point *x* in state (*p*, *t*) along the genome. The distribution of the distance to the first recombination encountered upstream (or downstream) from *x* is unaffected by drift. Thus the relationship between transition rates *Q* and discrete transition probabilities *R* is maintained: *Q*(*t*, *p* → *t*′, *p*′) = (*t* − 1)*R*(*t*, *p* → *t*′, *p*′) for (*t*, *p*) ≠ (*t*′, *p*′). If we indicate the state to the left or right of a recombination by a left- and right-pointing arrow, respectively, we write*c _{p}*

_{,}

_{t}_{,}

_{p}_{′}

_{,}

_{t}_{′}is the proportion of nonmigrant (diploid) sites with joint ancestry (

*p*,

*t*) and (

*p*′,

*t*′), which can be obtained using a recursive equation, as in Equation A3. In the driftless case, this reduces to

A case of particular interest is the pulse migration, with proportions *m* and 1 − *m* for populations 1 and 2, respectively. We then get *α*_{1}(*τ*) = *m*, and*N*→∞ yields the driftless case*T*→∞ reveals a linear dependence of the transition rate on the population size:

## Appendix 2: Numerical Estimation of Tract Length Distribution

In this section we describe how to obtain the expected distribution of tract lengths, given a set of Markov transition rates. A straightforward numerical solution strategy is to uniformize the transition matrix (Stewart 1994). Uniformization uses the fact that self-transition probabilities can be adjusted without affecting the trajectory statistics and in such a way that the total transition rate from each state is equal to the rate of the state with the highest transition rate, *Q*_{0}. Once all states have the same outgoing rate *Q*_{0}, the problem can be decomposed in two steps: a discrete calculation of the number of transitions in a given excursion and a calculation of the trajectory lengths given the number of transitions.

In the first step, we establish the distribution {*b _{n}*}

_{n}_{=1, …, ∞}of the number of steps spent in tracts of a given ancestry

*p*, which is a standard discrete Markov excursion problem. In principle, the number of steps can be arbitrarily large, but the probability of very long tracts decays rapidly, and after a certain number of steps the expected length of the excursion is more than the chromosome length. We therefore calculate {

*b*}

_{n}

_{n}_{=1,…,Λ}up to a cutoff Λ, such that

*Q*

_{0}>

*L*, the length of a chromosome). To ensure a proper probability distribution, we then set

*b*. For our purposes, we have found it convenient to evolve the state vector by repeated multiplication with a transition matrix modified to have a single, absorbing state corresponding to the non-

_{n}*p*ancestries and recording the amount of absorbed probability per multiplication.

The second step is straightforward since the length of the trajectories with *k* steps follows the Erlang distributionleading to the following expression for the tract-length distribution:

Ancestry tract-length distributions obtained in the infinite-chromosome limit may not be appropriate for finite genomes, particularly if many tracts have a length comparable to the chromosome length. For example, predicted tracts may be longer than the full chromosome length *L*, and these will not be observed. To model the tract-length distribution on a finite chromosome, we consider a general tract length distribution *φ*(*x*) on an infinite chromosome and ask for the distribution of tract lengths observed in a given window of length *L*. To this end, we first calculate the probability that the intersection of a tract of length *x*_{0} and a window of length *L* has length *x*. The probability *P*(*I*) that a tract of length *x*_{0} intersects the window of length *L* is proportional to *x*_{0} + *L*. Given *I*, and assuming that *x*_{0} < *L*, the probability that the intersection is of length *x* is*δ* Dirac’s delta function.

The result for *x*_{0} > *L* can be obtained by the permutation *x*_{0} ↔ *L*, so that*η* = min(*x*_{0}, *L*). This yields*L*, as*L* goes to infinity, the second term dominates and is proportional to *φ*(*x*).

The normalizing factor is*x* < *L* are considered. Finally, if only tracts of length >*C* are considered, we have

We now apply these results to the tract-length distributions from Equation A6. First, we note the Erlang distribution is related to the generalized incomplete gamma function by

## Appendix 3: Ancestry Variance in the Absence of Drift

### Ancestry variance under a Markov model of ancestry

We consider in this section the assortment variance, in the absence of drift, where ancestry in two individuals is modeled as independent realizations of a two-state Markov process. Let the Markov states representing ancestry be labeled by *k* = 1, 2 with rates *q*_{1} and *q*_{2} out of states 1 and 2, respectively. The generalization to multistate Markov processes is discussed below. We first consider a single chromosome of length *L* and are interested in the variance in *X*, the length of this chromosome covered in state *k* = 1. We have *ψ _{k}*(

*x*) the indicator function of state

*k*at position

*x*along the genome. Changing the order of the expectation and the integrals, we haveAll these expectations are independent of the position along the chromosome. We therefore havewith

*α*=

_{k}*q*

_{1−}

*/(*

_{k}*q*

_{1}+

*q*

_{2}) and

*P*(

_{k}*y*|

*x*) is the probability that

*y*is in ancestry

*k*given that

*x*is in ancestry

*k*. In a Markov process,

*q*

_{1}+

*q*

_{2}=

*T*− 2, and in the limit (

*q*

_{1}+

*q*

_{2})

*L*≫ 0, we recover our estimate

### Assortment variance for nonconstant migration

The generalization to arbitrary one-way migrations is straightforward, in the absence of drift. We evaluate Equation A16 by expanding on arrival times *s* for ancestry *p*:The probability that *y* is in the ancestry *p*, given that *x* is in state (*p*, *s*), can be written as *r* is the distance between *x* and *y* in morgans, ν represents a Markov state (*p*′, *s*′), *ρ _{νp}* is the indicator that

*p*′ =

*p*, and (

*κ*, υ

_{i}*) are the eigenvalues and eigenvectors of the transition matrix*

_{iν}*Q*

_{νν}_{′}. To obtain the

*s*in Equation A17 and solve only once for

*i*

_{0}with

### Distinguishing the two components of the ancestry variance from interchromosomal variance

We argued that, due to random chromosome assortment, the variance in ancestry between chromosomes is informative of the assortment variance. If all chromosomes had the same length, we could expect that the assortment variance in ancestry proportion across individuals would be proportional to the variance across chromosomes and inversely proportional to the number of chromosomes per individual. However, the different chromosomes have different lengths, and to combine the information we need an idea of how ancestry variance depends on chromosome length. We assume that the assortment variance on ancestry is inversely proportional to the chromosome length in morgans; in effect, we suppose that the number of independent ancestry observations is proportional to the chromosome length. The proportionality factor *σ _{g}* depends on the pedigree, so that

## Footnotes

*Communicating editor: M. A. Beaumont*

- Received February 17, 2012.
- Accepted March 28, 2012.

- Copyright © 2012 by the Genetics Society of America