About this Journal Submit a Manuscript Table of Contents
Advances in Bioinformatics
Volume 2008 (2008), Article ID 257864, 9 pages
http://dx.doi.org/10.1155/2008/257864
Review Article

A Tutorial of the Poisson Random Field Model in Population Genetics

1Department of Genetics, School of Medicine, University of Pennsylvania, Philadelphia, PA 19104, USA
2Department of Computer and Information Sciences, School of Engineering and Applied Sciences, University of Pennsylvania, Philadelphia, PA 19104, USA

Received 6 February 2008; Accepted 15 May 2008

Academic Editor: Alexander Zelikovsky

Copyright © 2008 Praveen Sethupathy and Sridhar Hannenhalli. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Population genetics is the study of allele frequency changes driven by various evolutionary forces such as mutation, natural selection, and random genetic drift. Although natural selection is widely recognized as a bona-fide phenomenon, the extent to which it drives evolution continues to remain unclear and controversial. Various qualitative techniques, or so-called “tests of neutrality”, have been introduced to detect signatures of natural selection. A decade and a half ago, Stanley Sawyer and Daniel Hartl provided a mathematical framework, referred to as the Poisson random field (PRF), with which to determine quantitatively the intensity of selection on a particular gene or genomic region. The recent availability of large-scale genetic polymorphism data has sparked widespread interest in genome-wide investigations of natural selection. To that end, the original PRF model is of particular interest for geneticists and evolutionary genomicists. In this article, we will provide a tutorial of the mathematical derivation of the original Sawyer and Hartl PRF model.

1. Introduction

Selectionists and neutralists have fiercely debated, for the past five decades, the extent to which Darwinian selection has shaped molecular evolution. However, both camps do agree that Darwinian selection is a bona fide natural phenomenon. Therefore, various so-called “tests of neutrality” have been developed to detect natural selection on a particular gene or genomic location (for a review on this topic, see [1]). However, these tests are often qualitative and only provide the directionality of selection. A decade and a half ago, S. Sawyer and D. Hartl provided a mathematical framework with which to determine quantitatively the intensity of selection on a particular gene, which they applied to the Adh locus in the Drosophila genome [2]. This framework is referred to as the Poisson random field (PRF) model. They then further used this framework to analyze codon bias in enteric bacteria [3]. Owing to the recent availability of whole genome sequences and genome-wide human polymorphism data, it has become increasingly tractable to perform genome-wide scans for signatures of selection. The PRF model has been applied to estimate the intensity of selection on synonymous and nonsynonymous sites throughout mitochondrial and nuclear genomes of a variety of species, including human [412]. Very recently, due to the advent of high-throughput experimental and computational identification of genomic regulatory elements, there has been an interest to estimate the intensity of natural selection on regulatory mutations. Chen and Rajewsky [13] use the PRF, among other techniques, to provide evidence for purifying selection (even stronger than on nonsynonymous coding sites) on a class of regulatory sites known as microRNA target sites. Due to the potentially wide range of applications of, and opportunities for theoretical extensions to, the PRF model, it is an increasingly important mathematical framework for quantitative geneticists. In this article, we will provide a tutorial of the mathematical derivation of the basic PRF model that was originally developed in [2]. The tutorial will follow the outline provided below: (i)Wright-Fisher model,(ii)diffusion approximation to the Wright-Fisher model,(iii)derivation, via diffusion theory, of formulas describing evolutionary processes of interest,(iv)derivation of the PRF using the above-mentioned formulas. The first three items are discussed in [14], and the last point was originally presented in [2]. In this tutorial, we aim to provide an integrated and comprehensive presentation that is accessible to nonprofessionals or beginners in the field of population genetics. Since the primary purpose is to review mathematical derivations, familiarity with calculus and at least a cursory knowledge of genetics will be helpful for the reader.

2. The Wright-Fisher Model

The Wright-Fisher (WF) model describes the change in frequency of a single mutation (derived allele) in a population over time. The simplest version of the model makes the following assumptions: (1) nonoverlapping generations, (2) constant population size in each generation, and (3) random mating, and is described as follows.

Consider a population of N diploid individuals that has a single polymorphic site with two alleles, one ancestral and one derived. Under this model, the frequency of the derived allele in the current generation is a function of the selection pressure on this allele and the binomial sampling effect with probabilities proportional to the frequency of this allele in the previous generation. The probability, pij, that there are j genes of the derived allele present at generation G + 1 given i genes of the derived allele present at generation G is given by the following binomial calculation: where depends on the relative fitness of the derived allele.

Assuming no dominance and no recurrent mutation, where is the fitness of the derived allele relative to 1 for the ancestral allele, and x (which is simply ) is the derived allele frequency (daf) in generation G. In the simplest model (no selection and no recurrent mutation), is simply x or .

The intuition behind Ψi is the following. Consider the scenario where both the ancestral and the derived alleles are neutrally evolving (no or negligible selection pressure). In this case, the probability of sampling a gene of the derived allele from the population in generation G is simply the frequency of the derived allele in generation G, or x. This can be rewritten as . Now, suppose that the derived allele is under some selection, s, meaning that the fitness of the derived allele is relative to 1 for the ancestral allele. In this case, genes are sampled according to their relative fitnesses (as in the equation for Ψi above). Figure 1(a) provides a pictorial representation of the basic Wright-Fisher model.

fig1
Figure 1: Pictorial representation of the Wright-Fisher process and its diffusion approximation: (a) basic Wright-Fisher model assuming selection, but no dominance or recurrent mutation and (b) diffusion approximation to the basic Wright-Fisher model.

3. Diffusion Theory

We define as the probability that a polymorphic site has i genes of the derived allele at time t, given that it had k genes of the derived allele at time 0. satisfies the following: where is given in (1).

It is convenient to change notation and write as , so that the above becomes In this framework, it has been shown to be extremely difficult to explicitly derive formulas for several quantities of evolutionary interest. However, as the size of the population approaches infinity (i.e., ), and assuming that the scaled selection pressure (Ns) and scaled mutation rate (Nμ) remain constant, the discrete Markov process given above can be closely approximated by a continuous-time, continuous-space diffusion process (Figure 1(b)): where is the probability distribution of x at time t, x is the daf at time t, p is the daf at time 0, and is the daf change in time .

We can perform a Taylor series expansion on both sides in and to derive the forward Kolmogorov equation: where and a(x) and b(x) depend on the genetic model (e.g., see eq (24)).

Equation (5) can be represented diagrammatically as in Figure 2. The probability of derived allele frequency at time is the product of the probability of moving from p to x in time t and the probability of moving from x to in time , summed over all possible values of x.

257864.fig.002
Figure 2: A diagrammatic intuition for (3) illustrates that the probability of derived allele frequency at time is the product of the probability of moving from p to x in time t and the probability of moving from x to in time , summed over all possible values of x.

The frequency trajectory of a derived allele can also be depicted as in Figure 3, which illustrates that the probability of frequency x at time is the product of the probability of moving from p to in time and the probability of moving from to x in time t, summed over all possible values of δp. This is formalized as follows: We can again perform a Taylor series expansion on both sides to derive the backward Kolmogorov equation: The forward and backward Kolmogorov equations have played a central role in theoretical population genetics since 1922. For details regarding their derivation, we refer the reader [15, Chapter 4]. Next, we will discuss how they are utilized to derive formulas for various quantities of evolutionary interest (yellow boxes in Figure 4).

257864.fig.003
Figure 3: A diagrammatic intuition for (5) illustrates that the probability of frequency x at time is the product of the probability of moving from p to in time and the probability of moving from to x in time t, summed over all possible values of .
257864.fig.004
Figure 4: The utility of the Kolmogorov equations for studying evolutionary processes. Blue boxes correspond to questions that clarify the assumptions of the genetic model being used, the red boxes correspond to when the Kolmogorov equations (KEs) are utilized, and yellow boxes correspond to quantities of evolutionary interest.

In a model where there is two-way recurrent mutation (i.e., there are no absorbing states, either extinction or fixation), stationarity is achieved when the probability of change in the derived allele frequency is no longer dependent on time t. We solve for the stationary distribution, f(x), in the following manner. First, we integrate through the forward Kolmogorov equation with respect to x : where is the probability of the derived allele assuming a frequency between 0 and x at time t. Therefore, the derivative of with respect to t can be interpreted as the probability flux (change in probability over time) of the diffusion process. The stationary distribution, f(x), can be solved by setting the probability flux equal to zero.

4. Derivation of Formulas Describing Evolutionary Processes of Interest

Let us now focus on a genetic model that assumes no recurrent mutation (i.e., two absorbing states, one at and another at ). As depicted by Figure 4, in such a model, it is possible to determine the probability of extinction (), the probability of fixation (), and the mean time until absorption (either at or ) by using the Kolmogorov backward equation (Figure 4). It is also possible to derive the mean time until absorption conditioned on always eventually reaching only one of the two states. Since this quantity is not directly applicable to the PRF, we do not review its derivation here, but instead refer the reader to [14].

4.1. Probability of Extinction

Using (11), we arrive at an equation parallel to (9): The probability that the derived allele frequency, x, reaches 0 at or before time t follows from (11) and is given by where p is the initial frequency of the derived allele and 0+ indicates 0 + ε, where ε is very small.

Replacing with , (12) can be written as As , can be interpreted as the probability that extinction ever occurs (independent of time) and can be rewritten in the form . From (14), it is evident that satisfies the following equation: Solving (15), we arrive at the following: where and where a(z) and b(z) are defined as in (6).

4.2. Probability of Fixation

The probability that the derived allele frequency, x, reaches 1 at time t follows from (11) and is given by where p is the initial frequency of the derived allele and 1 indicates , where ε is very small.

In (12), can be replaced by without any loss of generality. Also, by replacing with , (12) can be rewritten as By letting and solving for we arrive at the following: where has been defined in (17) and a(z) and b(z) have been defined in (6).

The probability of fixation and the probability of extinction must sum to 1. Using (16) and (20), we can verify that this is indeed the case.

Consider a genetic model that assumes the presence of selection, but no recurrent mutation, where and . Starting from (20), we can express the probability of fixation under this genetic model in the following manner:

4.3. Mean Time until Either Extinction or Fixation

We define to be the density function of the time t at which absorption occurs. The probability that absorption occurs, at either boundary or , by time t, is Furthermore, since absorption must happen by , we know that Performing integration by parts, we get the following: Equations (14), (19), and (22) show that satisfies the following equation: Using (25) and the fact that approaches 0 faster than t approaches , we can rewrite (24) as After interchanging the order of integration and differentiation we get where and is the mean time that the daf spends in the interval before absorption occurs.

We are interested in the case, where , since this is the initial frequency of the derived allele. In this case, we are interested only in values of x greater than , and for these values we can write and is defined in (17).

Under the simplest genetic model that assumes no selection and no recurrent mutation, we can set in (17) and (21) and show that reduces to p and reduces to 1. It follows from this that (29) can be reduced to Under a genetic model where , using , (29) can be rewritten as After integrating and simplifying the terms, we obtain Finally, substituting and , and invoking the approximation for small values of a, reduces approximately to where is a notation common in the literature to represent the expected time for which the population frequency of a derived allele is in the range before eventual absorption.

5. Poisson Random Field Theory

S. Sawyer and D. Hartl expanded the modeling of site evolution to multiple sites. Their model makes the following assumptions: (1) mutations arise at Poisson times, (2) each mutation occurs at a new site (infinite sites, irreversible), and (3) each mutant follows an independent WF process (no linkage). Sawyer and Hartl noticed from in (33), that is the expected number of sites in the population with derived allele frequency between x1 and (where equals , the per-locus mutation rate). The function , for which the full expression is given below, is also referred to in the literature as the limiting, equilibrium, or expected density function for derived allele frequencies. In a sample of size n, the expected number of sites with i (which ranges from 1 to ) copies of the derived allele is defined as a function of : The intuition behind F(i) is the following. The expected number of polymorphic sites with population daf x that have i copies of the derived allele out of n samples is given by the product of the expected number of sites with population daf , and the probability that each of those sites has i copies in the sample, which is given by the binomial calculation in the right-hand side of (36). To determine the expected number of sites with any population daf that have i copies of the derived allele, this product must be integrated over all possible values of x (resulting in F(i) above).

Consider the sample data , where is the observed number of sites with i copies of the derived allele out of n. Sawyer and Hartl showed that the number of derived alleles in the entire population at a particular frequency is a PRF with mean density given by (35) [2]. It follows, from the marking theorem on Poisson processes [16], that each random variable Xi is an independent Poisson distribution with mean equal to F(i) [2]. This framework allows us to define the probability of observing xi sites that have i copies of the derived allele (and copies of the ancestral allele) as the following: Since the Xi’s are independent, the probability of observing is given as The likelihood equation above provides a convenient means of estimating the values of the parameters θ and γ. The use of the PRF theory leads directly to a likelihood-ratio test of neutrality. Λ is defined as the ratio of the likelihood value under the maximum likelihood estimate of γ to the likelihood value under the neutral value of γ. It is a standard result that is asymptotically chi-square distributed with one degree of freedom [17].

Sawyer and Hartl further extended the PRF model in order to calculate the ratio of expected number of polymorphisms within species to expected number of fixed differences between species.In 1991, McDonald and Kreitman devised a 2-by-2 contingency table test of neutrality that was later named the MK test [18]. In the traditional MK test, a 2-by-2 contingency table is formed in order to compare the number of nonsynonymous and synonymous sites that are polymorphic within a species (RP and SP) and diverged between species (RF and SF) (Table 1). The central assumption of the MK test is that only nonsynonymous sites may be under selective pressure (i.e., synonymous sites are assumed to be neutrally evolving). If nonsynonymous sites are evolving according to a neutral model, then the expectation is that . However, if nonsynonymous sites are under negative selection, then the expectation is that , and if under positive selection, then . Sawyer and Hartl derived the formulas for the expected values of SP, SF, RP, and RF using their PRF theory [2]. Below are the derivations of each of these formulas. For all of the derivations, assume that the data consists of samples of size m and n from two different species.

tab1
Table 1: McDonald-Kreitman contingency table. 2-by-2 contingency table introduced by [18] for the inference of natural selection on nonsynonymous coding sites.
5.1. Expected Number of Synonymous Polymorphic Sites

Under neutral evolution (), the expected number of polymorphic sites with population daf x can be computed by taking the product of the per-locus mutation rate () and the probability under a neutral model of a single mutation having a frequency of x (from (30)): Now, consider species 1 with sample size m. The probability that a polymorphic site, with population daf equal to x, is detected as polymorphic in a sample of size m is given as The expected number of synonymous polymorphic sites, with population daf x, in the species 1 sample is the product of the expected number of synonymous polymorphic sites with daf x in the population () and the fraction of those that are expected to be detected in a sample of size . It follows then that the total expected number of synonymous polymorphic sites, with any population daf, in the species 1 sample is computed by integrating the product of and over the range of possible values for x: Finally, the total number of expected synonymous polymorphic sites in both species’ sample data is given as

5.2. Expected Number of Replacement Polymorphic Sites

The derivation of the expected value of RP follows the same logic. As described in (35), the expected number of polymorphic sites with population daf x given some average selection pressure γ is given by g(x). Similar to (41), the total expected number of replacement polymorphic sites in the species 1 sample is computed by integrating the product of g(x) and Pm(x) from 0 to 1: Finally, the total expected number of replacement polymorphic sites in both species’ sample data is given as

5.3. Expected Number of Synonymous Fixed Substitutions

When , the expected number of fixed substitutions in one species relative to another that diverged generations ago is given as the product of the number of total mutations and the probability of fixation of each mutation. The number of total mutations is the product of the mutation rate per generation and the number of generations since divergence is The probability of fixation is given in (21). As s approaches 0 (i.e., neutral evolution), the probability of fixationcan be reduced to p using the approximation for small values of a. Thus, for a newly derived neutral allele that has an initial frequency of , the probability of fixation is also .

Therefore, the total expected number of fixed substitutions in species 1 is However, given that the data are samples of the populations from both species, not all sites identified as fixed substitutions in the sample are truly fixed substitutions in the entire population. The expected number of sites in the species 1 sample that fall into this category is given by where and is given in (39).

Therefore, the total expected number of synonymous fixed substitutions in both species’ sample data is given as

5.4. Expected Number of Replacement Fixed Substitutions

Similar to the calculation of (46), given some selection pressure, γ, the expected number of fixed substitutions in one species relative to another that diverged generations ago is given as the product of (45) and (21): Substituting for p and invoking the approximation that for small values of a, we arrive at the following: However, again, given that the data are samples of the populations from both species, not all sites identified as fixed substitutions in the sample are truly fixed substitutions in the entire population. The expected number of sites in the species 1 sample that fall into this category is given by Therefore, the total expected number of replacement fixed substitutions in both species’ sample data is given as

5.5. Estimating Parameters

It is possible to obtain estimates of θ and γ by setting each of the observed values SP, RP, SF, and RF (Table 1) to their PRF expectations given by (42), (44), (48), and (52), respectively, and solving for the parameters. It has been shown that these estimates are equivalent to maximum-likelihood estimates [2, 19]. Bustamante et al. also eloquently describe and implement a hierarchical Bayesian model for parameter estimation [9].

6. Concluding Remarks

Sawyer and Hartl’s seminal presentation of the PRF in 1992 provided an innovative mathematical framework for estimating selection pressures and mutation rates, which are critical parameters that influence molecular evolution. However, it is worth noting that the model does harbor certain limitations. Foremost among these is the assumption of site independence, which is equivalent to the assumption of free recombination among mutations (i.e., no linkage). Thus, the model may not be appropriate for many data wherein strong linkage is present. Another limitation is the assumption of infinite sites (i.e., each mutation is at a new site). Although this assumption allows for a simpler model, it is not always biologically appropriate, especially for organisms that experience a higher mutation rate. Indeed, recent work has shown that the assumption of infinite sites can underestimate selection pressures and mutation rates and even infer positive selection, when in fact there is weak negative selection [20]. Recent theoretical work has focused on relaxing these and other assumptions of the original PRF model, so as to make it more appropriate for diverse biological contexts. For a brief list of such studies, we refer the reader to [20]. Ongoing theoretical and empirical work in this area will undoubtedly continue to extend the power of a PRF-based approach for population genetic inference.

Acknowledgments

The authors would like to thank Professors Joshua B. Plotkin and Warren J. Ewens for many thoughtful comments on the manuscript, discussions about the material, and suggestions about the presentation. Their support and expert advice have been instrumental to the successful completion of this tutorial. They also thank two anonymous reviewers for their helpful suggestions in improving the manuscript.

References

  1. S. Biswas and J. M. Akey, “Genomic insights into positive selection,” Trends in Genetics, vol. 22, no. 8, pp. 437–446, 2006. View at Publisher · View at Google Scholar · View at PubMed
  2. S. A. Sawyer and D. L. Hartl, “Population genetics of polymorphism and divergence,” Genetics, vol. 132, no. 4, pp. 1161–1176, 1992. View at Google Scholar
  3. D. L. Hartl, E. N. Moriyama, and S. A. Sawyer, “Selection intensity for codon bias,” Genetics, vol. 138, no. 1, pp. 227–234, 1994. View at Google Scholar
  4. H. Akashi, “Inferring weak selection from patterns of polymorphism and divergence at “silent” sites in Drosophila DNA,” Genetics, vol. 139, no. 2, pp. 1067–1076, 1995. View at Google Scholar
  5. M. W. Nachman, “Deleterious mutations in animal mitochondrial DNA,” Genetica, vol. 102-103, pp. 61–69, 1998. View at Publisher · View at Google Scholar
  6. D. M. Rand and L. M. Kann, “Mutation and selection at silent and replacement sites in the evolution of animal mitochondrial DNA,” Genetica, vol. 102-103, pp. 393–407, 1998. View at Publisher · View at Google Scholar
  7. H. Akashi, “Inferring the fitness effects of DNA mutations from polymorphism and divergence data: statistical power to detect directional selection under stationarity and free recombination,” Genetics, vol. 151, no. 1, pp. 221–238, 1999. View at Google Scholar
  8. D. M. Weinreich and D. M. Rand, “Contrasting patterns of nonneutral evolution in proteins encoded in nuclear and mitochondrial genomes,” Genetics, vol. 156, no. 1, pp. 385–399, 2000. View at Google Scholar
  9. C. D. Bustamante, R. Nielsen, S. A. Sawyer, K. M. Olsen, M. D. Purugganan, and D. L. Hartl, “The cost of inbreeding in Arabidopsis,” Nature, vol. 416, no. 6880, pp. 531–534, 2002. View at Publisher · View at Google Scholar · View at PubMed
  10. S. A. Sawyer, R. J. Kulathinal, C. D. Bustamante, and D. L. Hartl, “Bayesian analysis suggests that most amino acid replacements in Drosophila are driven by positive selection,” Journal of Molecular Evolution, vol. 57, pp. S154–S164, 2003. View at Publisher · View at Google Scholar · View at PubMed
  11. C. Bartolomé, X. Maside, S. Yi, A. L. Grant, and B. Charlesworth, “Patterns of selection on synonymous and nonsynonymous variants in Drosophila miranda,” Genetics, vol. 169, no. 3, pp. 1495–1507, 2005. View at Publisher · View at Google Scholar · View at PubMed
  12. C. D. Bustamante, A. Fledel-Alon, S. Williamson, et al., “Natural selection on protein-coding genes in the human genome,” Nature, vol. 437, no. 7062, pp. 1153–1157, 2005. View at Publisher · View at Google Scholar · View at PubMed
  13. K. Chen and N. Rajewsky, “Natural selection on human microRNA binding sites inferred from SNP data,” Nature Genetics, vol. 38, no. 12, pp. 1452–1456, 2006. View at Publisher · View at Google Scholar · View at PubMed
  14. W. J. Ewens, Mathematical Population Genetics: I. Theoretical Introduction, Springer, New York, NY, USA, 2004.
  15. W. J. Ewens, Mathematical Population Genetics, Springer, New York, NY, USA, 1979.
  16. J. F. C. Kingman, Poisson Processes, Oxford University Press, Oxford, UK, 1993.
  17. S. S. Wilks, Mathematical Statistics, John Wiley & Sons, New York, NY, USA, 1962.
  18. J. H. McDonald and M. Kreitman, “Adaptive protein evolution at the Adh locus in Drosophila,” Nature, vol. 351, no. 6328, pp. 652–654, 1991. View at Publisher · View at Google Scholar · View at PubMed
  19. S. Williamson, A. Fledel-Alon, and C. D. Bustamante, “Population genetics of polymorphism and divergence for diploid selection models with arbitrary dominance,” Genetics, vol. 168, no. 1, pp. 463–475, 2004. View at Publisher · View at Google Scholar · View at PubMed
  20. M. Desai and J. B. Plotkin, “Detecting directional selection from the polymorphism frequency spectrum,” Genetics, In press.