Monte Carlo methods for estimating population genetic parameters

1. Monte Carlo methods for estimating population genetic parameters Rasmus Nielsen University of Copenhagen

2. Outline • Idiosyncratic history and background on ML estimation of demographic parameters based on DNA sequence data. • A new computational approach/modification. • Idiosyncratic history and background on ML estimation of demographic parameters based on SNP data. • Ascertainment and large scale SNP data sets.

3. Felsenstein’s Equation So where Gi,i=1,2,…k, has been simulated from p(G|Q).

4. Coefficient of Variation

5. Importance Sampling So where Gi,i=1,2,…k, has been simulated from h(G).

6. Griffiths and Tavare Recursion Simulate mutation (coalescent) from and correct using importance sampling.

7. Example (Nielsen 1998) • Infinite sites model • Estimation of T • Estimation of population phylogenies

8. Integro-recursion Ugliest equation ever published in a biological journal…

9. MLE: T=1.8 (36,000 years)

10. Data from the Caribean Hawksbill Turtle

11. MCMC Set up a Markov chain on state space on all supported values of Q and G and with stationary distribution p(Q, G | X). Now since this can easily be done using Metropolis-Hastings sampling, i.e. updates to Q and G are proposed from a proposal distribution q(Q , G→ Q’ , G’) and accepted with probability

17. Some problems… • Histogram estimator or other smoothing must be used. • Likelihood ratios hard to estimate (e.g. M=0).

18. A new method • It is possible to calculate the marginal prior probability of a genealogy • It turns out that this math is doable, for most components of Θ such as q and M. • The we can sample from the marginal posterior of G using the previously discussed MCMC procedures. Slide inspired by Jody Hey

19. We then recover the posterior for Q using Approximated by Slide inspired by Jody Hey

20. Advantages • Eliminates problems with covariance between parameters leading to mixing problems. • Provides a smooth posterior/likelihood function useful for optimization and likelihood ratio estimation. Disadvantages • Requires more calculation in each MCMC iteration

21. Likelihood ratio estimation 6 loci, 15 gene copies, H0: m1=m2

22. Other approaches • Kuhner and Felsenstein use a combination of MCMC and importance sampling to estimate surfaces (no prior for the parameters). • PAC methods suggested by Stephens and Donnelly samples from a close approximation to to estimate an approximate likelihood. • ABC (Beaumont, Pritchard, Tavare and others) methods are a very popular and promising class of methods based on (1) reducing the data to summary statistics, (2) simulate new data from the prior, (3) accepting the parameter value under which the data was simulated if the difference between simulated and true statistics is less thand.

23. SNP DataNielsen and Slatkin (2000)

26. A more efficient method..Griffiths and Tavare (1998), Nielsen (2000)

27. A more efficient method..Griffiths and Tavare (1998), Nielsen (2000)

29. Ascertainment Sample vs. Typed Sample Ascertainment sample Typed sample

31. 1.0 0.9 0.8 E[D'] 0.7 0.6 0.5 0 1 2 3 4 5 6 7 8 9 10 r =2Nc no ascertainment bias ascertainment bias

32. Correcting for ascertainment biases Now, for simplicity, consider the case without a sweep, then where (in the simplest possible case) and

33. In this simple case, the maximum likelihood estimate of P is simply given by , k = 1, 2, …, n – 1, where nk is the number of SNPs with allele frequency k. Selective sweeps: Similarly define

34. 10,000 simulated SNPs with n = 20 and d = 5

35. Hudson’s (2001) Estimator when n = 100, m = 5, r = 5, and #SNP pairs = 200. Uncorrected Corrected

36. Complications • Double-hit ascertainment (HapMap) • Ascertainment based on chimpanzee (HapMap) • Panel depth may vary among SNPs and/or among regions (HapMap). • Ascertainment method may vary among SNPs (HapMap). • Population structure (HapMap). • Loss of information regarding asc. scheme (HapMap??).

37. HapMap ascertainment depth distrb. (ignores many important components)

39. Perlegen HapMap

40. Data Directly sequenced polymorphism data from 20 European-Americans, 19 African-Americans and one chimpanzee from 9,316 protein coding genes Data set previously described in Bustamante, C.D. et al. 2005. Natural selection on protein-coding genes in the human genome. Nature 437, 1153-7.

41. Demographic model European-Americans African-Americans migration Admixture Population growth Bottleneck

42. Estimation Sampling probabilities from the 2D frequency spectrum , Number of SNPs with pattern j in the 2D frequency spectrum SNPs within a gene are correlated. But estimator is consistent. The estimate has the same properties as a real likelihood estimator except that it converges slightly slower because of the correlation (Wiuf 2006).

43. African-Americans

44. European-Americans Godness-of-fit: p = 0.6

45. Acknowledgements Jody Hey, John Wakeley, Melissa Hubisz, Andy Clark, Carlos Bustamante, Scott Williamson, Aida Andres, Amit Andip, Adam Boyko, Anders Albrechtsen,Mark Adams, Michelle Cargill and other staff at Celera Genomics and Applied Biosystems.