Huff et al. BMC Genomics 2010, 11:8 http://www.biomedcentral.com/1471-2164/11/8 becomes fixed, or until t generations have passed. Let p equal the frequency of the allele in the subpopulation in the final generation. Then the importance weight for our desired distribution is the binomial likelihood func- tion:  n  x p(n  x)(1p). w    x  Because Wright-Fisher drift is a Markov process, the importance weight depends only on the allele frequency in the final generation. In contrast, Slatkin’s method employs a backward process that is only a rough approximation to Wright-Fisher drift, so the sampling weight must be calculated over the entire history of the two alleles with a separate term for each population and for each potential migration path in each generation [19]. Because the allele trajectory is generated from Wright- Fisher forward simulation, this method can seamlessly model complex demographic scenarios that include bot- tlenecks, expansions, and population subdivision with migration. The biggest downside to this flexibility is the potential for choosing parameter values that rarely result in population allele frequencies that are near the observed frequency in the sample. This concern must be evaluated when choosing parameter values, as some will require a prohibitive number of forward simulations to cover the sample space. However, all of the backward time methods are approximations to a forward Wright- Fisher process, and are meant to model natural pro- cesses that clearly occur in forward time, so this method is adequate for exploring most relevant models of posi- tive selection and demography. For models where the sample allele frequency is particularly unlikely, Slatkin’s method will be preferable since it involves a backward process conditioned on the sample [19]. For the results presented here, we set s and t to fixed values, though in principle they could be set to random variates in each forward iteration, reflecting uncertainty around estimates of selection strength and allele age. If t is a random variable, each origin generation-subpopula- tion must be weighted by its respective population size to reflect the probability that a new mutation originates in that generation [19]. Coalescent simulations We assumed all recombination events were crossovers, where a crossover occurs with the favored or neutral allele with probability proportional to the frequency of the alleles in the subpopulation [20,21]. For models with variable recombination rates, we followed Przeworski et al. [13] in setting the recombination rate to an exponen- tial random variate in each simulation, with mean equal to the rate of recombination in our constant models, 1 Page 7 of 9 cm/Mb. We adopted the implementation details of the coalescent process from [17], storing each generation in a lookup table indexed by the cumulative hazard of coa- lescence. To account for population subdivision, we introduced a subpopulation dimension to the coales- cence table. The trajectory of the favored allele was generated under a model where the migration rates are constant between subpopulations for each epoch. However, since a trajectory is in part a realization of this random pro- cess, we could not assume constant migration rates in a coalescent simulation based on a particular trajectory. The number of individuals of each genotype migrating to and from each population in a given generation is determined by the forward simulation and is therefore treated as a constant during the backward simulation. The individual migrants themselves are, however, cho- sen at random during the backward simulation. To implement this process, we introduced two migration lookup tables. The first table was analogous to the coa- lescence lookup table, storing the cumulative hazard of migration out of a given subpopulation for each allele. We used the second table to determine the destination subpopulation of a migrant, by storing the conditional probability of migrating from an origin subpopulation to a destination subpopulation given that a migration event occurred out of the origin subpopulation in a particular generation. Expanding on Coop and Griffith’s method, we accessed the coalescence and migration lookup tables with uniform random variates to generate the waiting time until the next event for each subpopula- tion-allele combination [17]. We then generate the wait- ing time until the next recombination event from an exponential random variate. Then from the memoryless property of the exponential distribution, the next event to occur is the event with the shortest wait time [17,20]. Ascertainment bias To introduce ascertainment bias to the simulated data, we developed a procedure to model the process in the Perlegen dataset. In their SNP discovery process, they identified all polymorphic sites in a fully sequenced sub- sample, then genotyped those sites in a larger sample [22]. To replicate their procedure on simulated data, we randomly assigned mutation events to the tree under an infinite sites model using a mutation rate of 2.2E-8 per nucleotide per generation [23]. We then designated 13- 33% of the chromosomes in each simulation as the ascer- tainment subsample [22]. We excluded all mutations that were not polymorphic in the subsample. To generate diploid genotypes for calculating FRC, we grouped the simulated chromosomes into randomly chosen pairs. Statistics EHH is defined as the probability that two chromo- somes in a sample share the same haplotype for a given