1. Introduction
Neutral loci linked to an allele of interest provide information about the history of that allele since it arose by mutation. The extent of linkage disequilibrium (LD) with linked markers allows tests of neutrality (Sabeti et al., Reference Sabeti, Reich, Higgins, Levine and Richter2002; Slatkin, Reference Slatkin2000; Slatkin & Bertorelle, Reference Slatkin and Bertorelle2001; Toomajian et al., Reference Toomajian, Ajioka, Jorde, Kushner and Kreitman2003; Voight et al., Reference Voight, Kudaravalli, Wen and Pritchard2006), estimates of selection intensity (Charlesworth, Reference Charlesworth2006; Charlesworth et al., Reference Charlesworth, Charlesworth and McVean2001; Slatkin, Reference Slatkin2001; Wiuf, Reference Wiuf2001a) and estimates of allele age (Goldstein et al., Reference Goldstein, Reich, Bradman, Usher and Seligsohn1999; Guo & Xiong, Reference Guo and Xiong1997; Kaplan et al., Reference Kaplan, Lewis and Weir1994; McPeek & Strahs, Reference McPeek and Strahs1999; Reich & Goldstein, Reference Reich, Goldstein, Goldstein and Schlötterer1999; Risch et al., Reference Risch, De Leon, Ozelius, Kramer and Almasy1995; Slatkin & Rannala, Reference Slatkin and Rannala1997; Stephens et al., Reference Stephens, Reich, Goldstein, Shin and Smith1998; Tishkoff et al., Reference Tishkoff, Dietzsch, Speed, Pakstis and Kidd1996). These problems are usually considered separately. When the goal is to estimate both age and selection intensity, age is estimated first and then selection intensity is inferred from the age by using a deterministic model of selection, an approach first taken by Stephens et al. (Reference Stephens, Reich, Goldstein, Shin and Smith1998).
In this paper, I consider the problem of jointly estimating the age of an allele and the intensity of selection it has experienced since it arose by mutation. I will formulate the problem in a Bayesian framework that accounts for the stochasticity of both allele frequency change and recombination. I then illustrate the approach by presenting a method for analysing low-frequency alleles linked to one or more marker loci. The method presented here is similar to that of Rannala & Reeve (Reference Rannala and Reeve2001) for linkage disequilibrium mapping in that it assumes that the dynamics of a low-frequency allele can be modelled by a linear birth–death process, and it approximates the probability of the data using Monte Carlo summation. It differs from Rannala & Reeve's method in using importance sampling instead of the Metropolis–Hasting algorithm for approximating the probability of the data.
One of the goals of this paper is to show that even under ideal conditions, when the genetic and demographic parameters are known without error, precise estimates of both allele age and selection intensity cannot be obtained because of the intrinsic unpredictability of genetic drift and recombination. This conclusion is in contrast to the conclusion reached when allele age is estimated using a ‘model-free’ method that assumes that the gene genealogy is a star (meaning that there are no internal branches) and that deterministic theory can be used to estimate the selection intensity under the assumption that the estimated age is the true age. The model-free approach does not take account of all sources of uncertainty in the estimates of allele age and selection intensity and hence exaggerates the accuracy of the resulting estimates.
2. Theory
(i) Formulation as a problem in Bayesian inference
Throughout, I will be concerned with the history of an allele, denoted by M, that arose as a unique mutation at time t 1 in the past in a population previously fixed for the alternative allele, m, at the locus of interest. All copies of M are identical by descent with the initial copy. Although t 1 is the true allele age, the time of the most recent common ancestor of all copies of M in a sample, called t 2 below, is also referred to as the age (Slatkin & Rannala, Reference Slatkin and Rannala2000). The goal is to estimate t 1 and to infer something about selection experienced by M since t 1. The genetic information available is the frequency of M in the population (x) and data from marker loci closely linked to M.
To describe the method in as simple a framework as possible, selection on M will be additive with selection coefficient s: the relative fitnesses of MM, Mm and mm individuals are 1+2s, 1+s and 1. I assume a prior distribution of s, Pr(s) that reflects what is known about selection on individual alleles. The strongest positive selection known in humans is on the S allele of the β-globin gene which causes sickle-cell anemia. Heterozygous carriers of S have roughly a 20% higher rate of surviving to adulthood than normal homozygous individuals in regions with a high incidence of malaria (Vogel & Motulsky, Reference Vogel and Motulsky1996). For analysing data from human populations, a reasonable prior distribution of s decreases with s and is small for s>0·2. In the analysis of G6PD, I assume Pr(s) = αe −αs, where α is a parameter with a value in the range 10–50.
The allele frequency, along with s and assumptions about past population growth, provide a prior distribution of t 1, Pr(t 1 |x, r, s) (Slatkin, Reference Slatkin2002b). Genetic data from linked marker loci provide additional information about s and t 1 that leads to the posterior distribution of s and t 1:
where G represents the genetic data.
How the posterior distribution is used depends on the goal of a study. If the goal is to estimate selection intensity, the best estimate of s is obtained by averaging over t 1, 
. If, instead, the goal is to estimate t 1 or to compare t 1 with historical information, then the posterior distribution of t 1 is obtained by averaging over s, 
. The resulting probability distribution can be used either to test a specific historical hypothesis, for example the probability that t 1 is larger or smaller than a given time, or to find the most probable value of t 1 and the associated confidence interval. In other situations, the goal is to find how much information about both t 1 and s is provided by the data, in which case the joint posterior distribution itself should be used.
(ii) Linear birth–death approximation for low-frequency alleles
To implement the above theory, I assume that the population has grown exponentially at rate r to its current size, N 0: the population size at time t in the past is N(t)=N 0e −rt. When x is small, the evolution of the number of copies of M can be approximated by a linear birth–death process (Wiuf, Reference Wiuf2001b), for which analytical expressions for Pr(t 1 |x, r, s) are available (Slatkin, Reference Slatkin2002b). The formulas used here are given in Appendix A.
The probability of the data given t 1, Pr(G|t 1, x, r, s), depends on the intra-allelic genealogy, meaning the gene tree of the i copies of M in a sample. As illustrated in Fig. 1, the intra-allelic genealogy is characterized by a set of i – 1 coalescence times, t 2, … , t i, and the topology, B (for branching pattern). The dependence of Pr(G| t 1x, r, s) on the coalescence times and topology can be written
where the sum is over all topologies with i tips and the integral is over all sets of intra-allelic coalescence times such that t 2 ⩾t 3 … ⩾t i >0 (Felsenstein, Reference Felsenstein1988). The first term in the integrand represents the effect of recombination on loci linked to M; it depends only on the intra-allelic genealogy and not on the probability attached to that genealogy. The second term represents the dependence of the joint distribution of intra-allelic coalescence times, which is easily derived for a linear birth–death process (Appendix A). Under the assumption that all copies of M are equivalent, the distribution of topologies is generated by assuming each lineage is equally likely to branch at each coalescent event, and hence depends only on i.
Fig. 1. The intra-allelic genealogy of a mutant M found in four copies (i=4), illustrating the definitions of allele age, t 1, and the intra-allelic coalescence times, t 2, t 3 and t 4.
(iii) Genetic data
I will consider here two kinds of genetic data. The first is the set of multi-locus haplotypes linked to each copy of M in the sample. I assume there are K diallelic marker loci, numbered 1 to K, with alleles at each locus labelled 0 and 1. The 1 allele is assumed ancestral in the sense that it was on the chromosome carrying the first copy of M. There are 2K possible haplotypes and their frequencies on non-M chromosomes are y j, j=1 … 2K, which are parameters of the model assumed to be known without error and not to have changed since t 1. The genetic data, G, consists of the haplotype associated with each M-bearing chromosome in the sample. I will discuss later the effects of uncertainly in y j and the haplotype phase of the markers on the M-bearing chromosomes.
The second kind of genetic data assumes that marker loci linked to M are sufficiently dense that it is possible to determine the length of the ancestral segment still linked to each copy of M. The frequency of each marker allele is unimportant. In this case the data consist of the lengths of ancestral segments measured in terms of the total recombination rate. I will assume the same rate c between every pair of sites. In that case, the data is a list of lengths in bases, l 1, … , l i, of ancestral segments linked to each copy of M in the sample. This kind of data differs from the minimum length of the ancestral fragments linked to M, l=min[l 1, … , l i], which has been used as a summary statistic (McPeek & Strahs, Reference McPeek and Strahs1999; Slatkin, Reference Slatkin2001; Slatkin & Bertorelle, Reference Slatkin and Bertorelle2001).
For both kinds of data, the computational problem is to find Pr(G| t 1, x, r, s). As is usually the case when summing over tree topologies and integrating over branch lengths, exact computations are not feasible for more than three or four M-bearing chromosome. Here, I use an importance sampling method introduced previously (Slatkin, Reference Slatkin2002a). The transition matrices used for the two types of data are summarized in Appendix B. The results from using the importance sampling method were checked by comparing them with results for one and two marker loci which could be obtained by other means.
3. Positive or negative selection
The Bayesian framework described in the previous sections does not require that selection favoured of the mutant of interest. In fact, unless an additional step is taken, it is not possible to determine whether the mutant experienced positive or negative selection. The problem is seen most clearly in the case of an allele with an additive effect on fitness in a population of constant size. Maruyama (Reference Maruyama1974) showed that, under these assumptions, the prior distribution of age is independent of the sign of the selection coefficient. Although Maruyama's conclusion is no longer true if there is dominance or past population growth, the implication of his result for the problem of estimating the selection coefficient from allele age is still important. A deleterious allele, if it is found in any substantial frequency, has to be relatively young and hence in LD with linked markers, just as would an advantageous allele.
Maruyama's result led Nordborg & Tavaré (Reference Nordborg and Tavaré2002) to conclude that extensive LD near an allele might well indicate selection against it rather than in favour of it. They suggested that which is more probable depends on the relative rates of deleterious and advantageous mutations. And they concluded that our current lack of knowledge about mutation parameters renders the problem of deciding whether selection was in favour of or against a mutant ‘philosophical’. Nordborg & Tavaré's (Reference Nordborg and Tavaré2002) argument is not convincing, however, because it ignores the fact that an advantageous mutation is much more likely to reach a given frequency than will a mutation that is deleterious to the same extent. In Appendix C, I show how to use the birth–death approximation to find the ratio of probabilities that a mutant was advantageous instead of deleterious. If the mutant is present in any substantial frequency, that ratio is so large that the possibility of negative selection can be ignored. For example, in the case of G6PD, discussed below, it is approximately 1·4×1016 more likely that the selection coefficient affecting the A− allele of G6PD is +0·05 than −0·05.
4. Application to G6PD
The Bayesian method is illustrated by applying it to the data of Saunders et al. (Reference Saunders, Slatkin, Garner, Hammer and Nachman2005) who sequenced 20 chromosomes carrying the A− allele at the X-linked G6PD locus. This allele is found in sub-Saharan African populations in frequencies up to 20%. It is thought to confer partial resistance to malaria (Vogel & Motulsky, Reference Vogel and Motulsky1996). Saunders et al. (Reference Saunders, Slatkin, Garner, Hammer and Nachman2005) sequenced eight loci in a 2·5 Mb region surrounding G6PD in 51 unrelated males of sub-Saharan African origin. Twenty of these males carried the A− allele. LD at linked SNP markers spanned a region of 1·6 Mb on the A− chromosomes. By applying a simplified version of the method presented here, Saunders et al. (Reference Saunders, Slatkin, Garner, Hammer and Nachman2005) concluded that heterozygous carriers of A− had between a 10% and 20% fitness advantage over homozygous individuals and that A− arose by mutation less than 150 generations (3750 years) ago. Part of the purpose of this paper is to formalize and generalize the method used in that paper and in the paper of Wood et al. (Reference Wood, Stover, Slatkin, Nachman and Hammer2005).
Following Saunders et al. (Reference Saunders, Slatkin, Garner, Hammer and Nachman2005), I assume no population growth (r=0) and a frequency, x=0·11, of A−. With a population size of N 0=100 000, a fraction f=9·1×10−4 of the population would have to be sampled to obtain 20 copies of A− in a sample.
To analyse the linked haplotypes, I considered four SNPs (numbers 36, 41, 55 and 99), one each from the four loci BGN, IDH3G, 1CAM and G0·9MT in figure 2 of Saunders et al. (Reference Saunders, Slatkin, Garner, Hammer and Nachman2005). I assumed that the nucleotide carried by the first individual listed (ALB77) was on the ancestral A− chromosome and denoted it by 1. The other nucleotide was denoted by 0. With this convention, each haplotype can be numbered from 0 to 15 by treating the haplotype as a binary number. With this convention, the data of Saunders et al. can be reduced to six A− chromosomes with haplotype 15, three with haplotype 9, seven with haplotype 7, one with haplotype 3 and three with haplotype 1. I used the haplotypes of the 31 non-A− chromosomes to estimate the frequencies of the 1 allele to be 0·419, 0·0323, 0·0645 and 0·433 for the four loci. There is no obvious LD on the non-A− chromosomes, so I computed the background frequencies of the 15 haplotypes by assuming linkage equilibrium.
The linkage map for these four loci was inferred by assuming that the rate per base is constant in this region and using the estimated rate between L1CAM and G6PD of 0·01675 to estimate the rate per base. The map distances from BGN were 0·0086 for IDH3G, 0·0131 for L1CAM, 0·02985 for G6PD and 0·0354 for G0·9MT. From the background frequencies and map distances, the elements of the 16×16 transition matrix were obtained.
The other way of analysing this kind of data is to infer the length of the ancestral chromosome still linked to A−. That assessment is somewhat subjective. For example, in the G6PD data set the first four of the A− chromosomes (ALB77, SHO07, VA065 and IVC17) have nearly the same haplotype as far as the BGN locus, which is 991 kb from G6PD. But the sequence of VA065 differs from the other three at SNP 20 and the sequence of IVC17 differs from the other three at SNP 36, both of which are in the BGN locus. The two apparently aberrant sites could differ from the others because of gene conversion or mutation that did not otherwise affect the ancestral fragment. Or they could indicate that recombination had taken place between BGN and IDH3G and then subsequent recombination recreated what appears to be the ancestral haplotype. I assumed that a single aberrant SNP was not sufficient evidence to indicate recombination. With that assumption, the data of Saunders et al. (Reference Saunders, Slatkin, Garner, Hammer and Nachman2005) reduce to the following: four of the A− chromosomes had the conserved ancestral fragment to BGN on the centromeric side of G6PB, seven to IDH3G, one to L1CAM, four to TAZ and one had no conserved fragment. I did not analyse the telomeric side of G6PD because all 20 A− chromosomes had identical haplotypes on that side. I rounded the distances to 50 bp intervals. BGN was assumed to be 20 units from G6PD, IDH3G was 14 units, L1CAM was 11 units and TAZ was 2 units. The recombination distance per 50 bases was 0·01675/11=0·001528, based on the estimated recombination rate between G6PD and L1CAM.
The results from analysing the two data sets are presented in Fig. 2. Part (a) shows the joint posterior distribution of s and t 1 based on the analysis of ancestral fragments; (b) shows the same distribution obtained from the analysis of four marker loci. The two distributions have the same character. There is a ridge indicating that t 1 and s are somewhat confounded. The reason is that both the prior distribution of allele age and the distribution of intra-allelic coalescence times depend strongly, although not exclusively. on the product st 1. The two data sets nevertheless result in joint estimates of s and t 1 that are consistent with each other. The maximum for Fig. 2a is at s=0·26 and t 1=40 generations and for Fig. 2b at s=0·24 and t 1=40 generations.
Fig. 2. The joint posterior distribution of s and t 1 computed for the A− allele of G6PD using the importance sampling method described in the text. The data and genetic map were taken from Saunders et al. (Reference Saunders, Slatkin, Garner, Hammer and Nachman2005). In both cases, 106 replicates were used. (a) Pr(s, t 1) estimated from the lengths of ancestral chromosomal segments linked to A−. The identification of the ancestral segments is described in the text. (b) Pr(s, t 1) estimated from data for four marker loci (BGN, IDH3G, 1CAM and G0·9MT) linked to G6PD. The details of the genetic map are described in the text.
The similarity of the most probable values of s and t 1 does not tell the whole story. Although it is somewhat difficult to see in the three-dimensional graphs in Fig. 2, there is much more variability in Fig 2b than in 2a. The variability is more obvious when the marginal distributions of s and t 1 are computed (Fig. 3). The marginal distributions of s and t 1 for the model of ancestral segments (Fig. 3a and b) and are much smoother than the corresponding distributions based on the four marker loci (Fig. 3c and d). The results for both models are based on 106 replicates of the importance sampling algorithm, so the difference between results from the two models reflects the greater intrinsic variability in recombination in a small set of marker loci.
Discussion
The Bayesian framework for jointly estimating allele age and selection intensity is quite general and can be adapted to other models of selection and to other assumptions about population history. The implementation of the method presented here is appropriate for low-frequency alleles. The application to the A− allele of G6PD shows that allele age and selection intensity are partially confounded because of the underlying population genetic processes. We can see why that is true by considering a deterministic model of selection on an allele with additive effect s on fitness. The differential equation that approximates the change in allele frequency is
where p is the allele frequency and t is time (Hartl & Clark, Reference Hartl and Clark1997). This equation is unchanged if s is multiplied by a constant C, sˈ = Cs, and the units of time are reduced by the same factor, tˈ = t/C. Hence any function that depends on the product st is unchanged.
The reason that s and t 1 are not completely confounded is that, when M is at low frequency, the trajectory of allele frequencies is not deterministic; it is also affected by genetic drift. Drift and selection are both accounted for in the linear birth–death approximation. The equations in Appendix A show that both the prior distribution of age and the joint distribution of intra-allelic coalescence times depend on s and t 1 separately as well as on the product st 1.
The Bayesian method accounts for both the stochastic dynamics of a low-frequency allele and the stochastic nature of recombination. These two sources of unpredictability make it difficult to obtain precise estimates of age and selection intensity even under the idealized conditions assumed here. Any additional sources of uncertainty, for example in estimated map distances, estimated haplotype frequencies of the non-M chromosomes or inferred haplotype phases, further reduce the ability to estimate ages and selection intensities.
The importance of accounting for all sources of variability is evident in the analysis of the G6PD data set. Figs 2 and 3 show that there is more variability in the results for the four marker loci than for the ancestral segments. The reason is that, in the model of individual marker loci, there are more sources of uncertainty. In this example, the age of the A− allele is small enough (≈50 generations) and the recombination distances between the markers are also small enough (≈1–3 cM) that very few recombination events are expected after A− arose by mutation. The data do not allow a precise determination of the number or order of those recombination events and hence do not narrowly constrain the estimates of s and t 1. Failure to take account of all sources of variability, as in ‘model-free’ methods that assume that the intra-allelic genealogy is a star (Neuhausen et al., Reference Neuhausen, Mazoyer, Friedman, Stratton and Offit1996; Reich & Goldstein, Reference Reich, Goldstein, Goldstein and Schlötterer1999), lead to exaggerated confidence in the resulting estimates.
The estimated selection intensity on the A− allele of G6PD, 0·24, is larger than the minimum estimate of 0·05 reported by Saunders et al. (Reference Saunders, Slatkin, Garner, Hammer and Nachman2005), and the estimated age, 40 generations, is smaller than their estimate. The difference is a consequence of assuming a prior distribution for s and using it in the Bayesian analysis. Saunders et al. (Reference Saunders, Slatkin, Garner, Hammer and Nachman2005) found that the likelihood function depends only weakly on s for s>0·05. They used s=0·05 as a conservative estimate and then estimated t 1 from the posterior distribution of t 1, given that value of s. The result is an overestimate of t 1 and an underestimate of s. In contrast, the estimates obtained here are from a posterior distribution that has a single maximum, albeit one that lies on an evident ridge.
In the G6PD example, no population growth was assumed. If there actually had been exponential growth, then the estimated selection coefficient would estimate r+s because population growth has almost the same effect on linked markers as does positive selection. However, exponential growth would affect all loci in the same way, resulting in more LD across the genome. The evidence for positive selection at G6PD comes from the fact that the extent of LD with the A− locus is unusually large.
In the general formulation of the birth–death approximation, only exponential growth at rate r was considered in order to order to minimize the number of parameters in the model. The method could be adapted to other models of demographic history, including population subdivision and bottlenecks in population size. As in the case of exponential growth, evidence of selection comes from LD with the allele of interest that exceeds background levels of LD created by demographic processes.
The Bayesian analysis described above assumes a prior distribution of selection intensities, but treats the observed mutant frequency (x) as a fixed parameter. It would be possible to develop a Bayesian method in which a prior distribution of x is derived from assumptions about mutation, selection, drift and other factors. Such an approach, however, would not take account of the way that mutants are chosen for detailed genetic analysis of the type done by Saunders et al. (Reference Saunders, Slatkin, Garner, Hammer and Nachman2005). Until the criteria for carrying out such studies can be quantified, it seems difficult to extend the Bayesian framework to include a prior distribution of x as well.
This research was supported in part by Grant GM40282 from the US National Institutes of Health. I thank B. Rannala for helpful discussions of this topic and the reviewers for helpful comments on a previous version of this paper.
Appendix A. Approximations based on the linear birth–death process
Prior distribution of allele age
The prior distribution is itself derived from a Bayesian formulation initially used by Kimura & Ohta (Reference Kimura and Ohta1973). The probability that an allele M of frequency x 0 at t=0 is at frequency x at a later time t=T can be written as a transition function ɸ(x, T | x 0, 0) that depends on the population sizes and selection intensities between 0 and T. The prior distribution of allele age is obtained from this transition function by assuming a constant mutation rate μ. The probability that a mutation arises t 1 generations in the past is 2 μN(−t 1) where N(−t 1) is the population size at that time. The probability that M arises t 1 generations in the past and is found at frequency x today is
where x 0 is the frequency of a new mutant, 1/[2N(−t 1)]. Therefore,
For neutral alleles in a population of variable size and for arbitrary x 0 and x, ɸ can be obtained as the solution to a diffusion equation that approximates the Wright–Fisher, Moran and other models of gene frequency evolution. For selected alleles, there is no usable solution for ɸ even in a population of constant size. If x 0 and x are both small, the problem can be simplified because the number of copies of M can be modelled by a linear birth–death process that allows for both exponential population growth and additive selection. (Slatkin, Reference Slatkin2002b; Wiuf, Reference Wiuf2001b). In this case, analytical expressions for Pr(t 1) are available. If there are N 0 individuals in the population at the present time and n chromosomes are chosen at random and tested for M, the fraction of the population sampled is f=n/(2N 0) and the number of copies of M in the sample is i=2N 0x. The population is assumed to have grown exponentially at rate r and s is the additive selection coefficient in favour of M. The growth rate r is non-negative, but s may take any value.
I have presented elsewhere distributions for all combinations of r and s (Slatkin, Reference Slatkin2002b). Here I include only those formulas used to obtain the results presented in the text. If r=s=0, the prior probability that M arose t generations in the past is
if r=0 and s>0,
if r>0 and s>0,
where ξ=r+s, Γ( ) is the gamma function and 2F 1( ) is the hypergeometric function (Abramowitz & Stegun, Reference Abramowitz and Stegun1965). In these equations, t is used instead of t 1 for notational convenience. Note that in equations (A2b) and (A2c), Pr(t) depends on s and t separately as well as on the product st.
Intra-allelic coalescence times
As shown by Slatkin & Rannala (Reference Slatkin and Rannala1997), the coalescence times for the intra-allelic genealogy are easy to generate if the number of copies of the mutant follows a linear birth–death process. A random set of coalescence times {t 2, … , t i} is obtained by generating i−1 random variables from a kernel distribution, b (t| t 1), and then sorting them in decreasing order. If r=s=0,
and if ξ=r+s>0,
If the intra-allelic genealogy were a star, then t 2=…=t i=t 1. In this model, that would be equivalent to b (t| t 1) being a spike (i. e. a Dirac δ function) at t=t 1. The extent to which b (t| t 1) is not a spike indicates how much the intra-allelic genealogy differs from a star.
Appendix B. Estimating the probability of the data by using importance sampling
The method described in the text requires the calculation of Pr(G| t 1, x, r, s), which is the probability of the genetic data (G), given the age of M (t 1), the frequency of M (x), the population growth rate (r) and the selection intensity (s). The exact calculation requires the evaluation of equation (2) in the text. Because the number of branching patterns of a genealogy increases very rapidly with the number of terminal branches, it becomes impractical to evaluate the sum exactly for a sample size larger than 3 or 4. To avoid this problem, I use an importance sampling method introduced previously (Slatkin, Reference Slatkin2002a) that performs well when the intra-allelic genealogy is generated by a linear birth–death process. The idea is to sample branching patterns in such a way that patterns that contribute more to the overall probability of the data are sampled more frequently. By generating a large number of trees and taking the average, an estimate of the probability of the data is obtained:
where H is the total number of replicates, PrRS(B h) is the probability of branching pattern B h under a random choice of branching patterns, PrIS(B h) is the probability of B h under the importance sampling method used and F is the forward transition matrix for the genetic data in one generation. This transition matrix has to be specified for each type of data. In this paper, I consider two models: one with K loci and an arbitrary recombination map and the other of the length of the ancestral chromosome associated with each copy of M.
Multiple marker loci
There are assumed to be K diallelic marker loci linked to M. At each locus, allele 0 is assumed to have been on the ancestral M chromosome. The 2K haplotype frequencies on the non-M chromosomes are assumed to be known and fixed. The genetic map is specified by c, the recombination rate between M and locus 1. If c<0, locus 1 is to the left of M; if c>0, locus 1 is the to right. The recombination rate between locus k and k+1 is c k (k=1, … , K−1). Therefore, M is the leftmost locus if c>0 and M is the rightmost locus if 
. Otherwise M is between one pair of marker loci. The number of M-bearing chromosomes is assumed to be so small that only recombination events between an M-bearing and a non-M-bearing chromosome are considered. There is no interference in recombination and the recombination rates are assumed to be sufficiently small that at most one crossing-over occurs each generation.
With these assumptions, it is straightforward to find the elements of F, the forward transition matrix for the 2K states of the marker loci on an M-bearing chromosomes. Each off-diagonal element of F represents a transition from one marker haplotype to another. The transition probability is the sum of the probabilities of recombination events between each pair of marker loci. For example, if K=4, the probability of a transition from (1, 0, 1, 0) to (1, 0, 1, 1) can result from a recombination event between markers 2 and 3 bringing in the 1, 1 haplotype at markers 3 and 4 and a recombination event between markers 3 and 4 bringing in the 1 allele at marker 4. The diagonal elements are obtained by subtracting from 1.
Ancestral chromosome lengths
If the lengths of ancestral chromosomes linked to each copy of M are known, it is possible to derive an analytical expression for an arbitrary power of F, the forward transition matrix. The chromosomal segments on each side of M are considered separately and then the probabilities are multiplied to obtain the overall probability. On one side of M, assume there are L sites and the recombination rate between each pair of sites is c. The parameter L is chosen to be large enough that the results do not depend on its value. The state of a segment is k, the length of the ancestral segment linked to M (k=0, 1, … , L).
The elements of F are probabilities of transitions from state k in generation t to state k′ in generation t+1: F kk′ = 0 if k′>k, F kk′ = (1 − c)k if k′=k and F kk′ = c (1 − c)k′ if k′<k. It is straightforward to verify that the elements of F t are 0 if k′>k, (1 − c)kt if k′=k and (1 − c)k′t − (1 − c)(k′ + 1)t if k′<k.
Appendix C. Positive or negative selection
The theory of linear birth–death processes provides the probability that j alleles are found at time t, given one copy at time 0, a birth rate λ and a death rate μ:
(j>0), where
(Kendall, Reference Kendall1948). If a fraction f of the population is sampled with replacement, then the probability that i copies of the allele are present in the sample is obtained by multiplying (C1) by a binomial distribution with probability f and sample size j and summing over j to obtain:
where u′ = fu/[1 − (1 − f)u].
Slatkin & Rannala (Reference Slatkin and Rannala1997) and Wiuf (Reference Wiuf2001b) showed that the linear birth–death process approximates the dynamics of a rare allele in a population of size N growing exponentially at rate r and subject to additive selection of intensity s if λ=1/2 and μ=1/2−r−s.
Given that i copies of a mutant are observed at time t, the relative likelihood that a mutant has a selective advantage s instead a selective disadvantage of s is the ratio Pr(i|s, t)/Pr(i| − s, t). This ratio is very large if s is not small. In the example used in the text, in which i=20, f=9·1×10−4 and r=0, the ratio is approximately 1·4×1016 if s=0·05 and increases rapidly as s increases.
