A Tutorial of the Poisson Random Field Model in Population Genetics

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 [4–12]. 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, p_ij, 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: 𝑝_𝑖𝑗=𝑗2𝑁(Ψ_𝑖)^𝑗(1−Ψ_𝑖)^2𝑁−𝑗, (1) where Ψ_𝑖 depends on the relative fitness of the derived allele.

Assuming no dominance and no recurrent mutation, Ψ_𝑖=𝑥(1+𝑠),𝑥(1+𝑠)+(1−𝑥) (2) where 1+𝑠 is the fitness of the derived allele relative to 1 for the ancestral allele, and x (which is simply 𝑖/2𝑁) is the derived allele frequency (daf) in generation G. In the simplest model (no selection and no recurrent mutation), Ψ_𝑖 is simply x or 𝑖/2𝑁.

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, 𝑖/2𝑁 or x. This can be rewritten as 𝑥/[𝑥+(1−𝑥)]. Now, suppose that the derived allele is under some selection, s, meaning that the fitness of the derived allele is 1+𝑠 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.

(a)

(b)

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: 𝑝^(𝑡+1)_𝑘𝑗=_𝑖𝑝^(𝑡)_𝑘𝑖𝑝_𝑖𝑗, (3) where 𝑝_𝑖𝑗 is given in (1).

It is convenient to change notation and write 𝑝^(𝑡)_𝑘𝑖 as 𝑓(𝑥;𝑝,𝑡), so that the above becomes 𝑓(𝑗;𝑘,𝑡+1)=_𝑖𝑓(𝑖;𝑘,𝑡)𝑝_𝑖𝑗. (4) 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)): 𝑓(𝑥+𝛿𝑥;𝑝,𝑡+𝛿𝑡)=¹₀𝑓(𝑦;𝑝,𝑡)𝑓(𝑥+𝛿𝑥;𝑦,𝛿𝑡)𝑑𝑦, (5) 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: 𝜕𝑓(𝑥;𝑝,𝑡)=𝜕𝜕𝑡²𝑏(𝑥)𝑓(𝑥;𝑝,𝑡)2𝜕𝑥²−𝜕𝑎(𝑥)𝑓(𝑥;𝑝,𝑡),𝜕𝑥 (6) where 𝐸(Δ𝑥)≈𝑎(𝑥)𝑑𝑡,var(Δ𝑥)≈𝑏(𝑥)𝑑𝑡, (7) 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.

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: 𝑓(𝑥;𝑝,𝑡+𝛿𝑡)=¹₀𝑓(𝑝+𝛿𝑝;𝑝,𝛿𝑡)𝑓(𝑥;𝑝+𝛿𝑝,𝑡)𝑑(𝛿𝑝). (8) We can again perform a Taylor series expansion on both sides to derive the backward Kolmogorov equation: 𝜕𝑓(𝑥;𝑝,𝑡)𝜕𝜕𝑡=𝑏(𝑝)²𝑓(𝑥;𝑝,𝑡)2𝜕𝑝²𝜕+𝑎(𝑝)𝑓(𝑥;𝑝,𝑡).𝜕𝑝 (9) 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).

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 𝛿𝑝.

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 : 𝜕𝐹(𝑥;𝑝,𝑡)=𝜕𝜕𝑡𝑏(𝑥)𝑓(𝑥;𝑝,𝑡)−,2𝜕𝑥𝑎(𝑥)𝑓(𝑥;𝑝,𝑡)𝐹(𝑥;𝑝,𝑡)=^𝑥₀𝑓(𝑦;𝑝,𝑡)𝑑𝑦, (10)(11) 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 𝑥=0 and another at 𝑥=1). As depicted by Figure 4, in such a model, it is possible to determine the probability of extinction (𝑥=0), the probability of fixation (𝑥=1), and the mean time until absorption (either at 𝑥=0 or 𝑥=1) 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): 𝜕𝐹(𝑥;𝑝,𝑡)𝜕𝜕𝑡=𝑏(𝑝)²𝐹(𝑥;𝑝,𝑡)2𝜕𝑝²𝜕+𝑎(𝑝)𝐹(𝑥;𝑝,𝑡).𝜕𝑝 (12) The probability that the derived allele frequency, x, reaches 0 at or before time t follows from (11) and is given by 𝑃₀(𝑝,𝑡)=⁰⁺₀𝑓(𝑦;𝑝,𝑡)𝑑𝑦=𝐹0⁺;𝑝,𝑡, (13) where p is the initial frequency of the derived allele and 0⁺ indicates 0 + ε, where ε is very small.

Replacing 𝐹(0⁺;𝑝,𝑡) with 𝑃₀(𝑝,𝑡), (12) can be written as 𝜕𝑃₀(𝑝,𝑡)𝜕𝜕𝑡=𝑏(𝑝)²𝑃₀(𝑝,𝑡)2𝜕𝑝²𝜕𝑃+𝑎(𝑝)⁰(𝑝,𝑡).𝜕𝑝 (14) 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: 𝜕0=𝑏(𝑝)²𝑃₀(𝑝)2𝜕𝑝²𝜕𝑃+𝑎(𝑝)⁰(𝑝).𝜕𝑝 (15) Solving (15), we arrive at the following: 𝑃₀∫(𝑝)=^1𝑝𝜓(𝑦)𝑑𝑦∫¹₀,𝜓(𝑦)𝑑𝑦 (16) where 𝜓(𝑦)=𝑒^{−2∫𝑦[𝑎(𝑧)/𝑏(𝑧)]𝑑𝑧} (17) 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 𝑃₁(𝑝,𝑡)=¹₁₋𝑓(𝑦;𝑝,𝑡)𝑑𝑦=1−¹⁻₀𝑓(𝑦;𝑝,𝑡)𝑑𝑦=1−𝐹(1⁻;𝑝,𝑡), (18) where p is the initial frequency of the derived allele and 1⁻ indicates 1−𝜀, where ε is very small.

In (12), 𝐹(𝑥;𝑝,𝑡) can be replaced by 1−𝐹(𝑥;𝑝,𝑡) without any loss of generality. Also, by replacing 1−𝐹(1⁻;𝑝,𝑡) with 𝑃₁(𝑝,𝑡), (12) can be rewritten as 𝜕𝑃₁(𝑝,𝑡)𝜕𝜕𝑡=𝑏(𝑝)²𝑃₁(𝑝,𝑡)2𝜕𝑝²𝜕𝑃+𝑎(𝑝)¹(𝑝,𝑡).𝜕𝑝 (19) By letting 𝑡→∞ and solving for 𝑃₁(𝑝), we arrive at the following: 𝑃₁∫(𝑝)=^𝑝0𝜓(𝑦)𝑑𝑦∫¹₀,𝜓(𝑦)𝑑𝑦 (20) 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 𝑎(𝑥)=𝑠𝑥(1−𝑥) and 𝑏(𝑥)=𝑥(1−𝑥)/2𝑁. Starting from (20), we can express the probability of fixation under this genetic model in the following manner: 𝑃₁∫(𝑝)=^𝑝0𝑒^{−2∫𝑦[𝑎(𝑧)/𝑏(𝑧)]𝑑𝑧}𝑑𝑦∫¹₀𝑒^{−2∫𝑦[𝑎(𝑧)/𝑏(𝑧)]𝑑𝑧}=∫𝑑𝑦_𝑝0𝑒^{−4𝑁𝑠𝑦}𝑑𝑦∫¹₀𝑒^{−4𝑁𝑠𝑦}=𝑑𝑦1−𝑒^{−4𝑁𝑠𝑝}1−𝑒^−4𝑁𝑠. (21)

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 𝑥=0 or 𝑥=1, by time t, is 𝑃₀(𝑝,𝑡)+𝑃₁(𝑝,𝑡)=^𝑡₀𝜙(𝑝,𝑡)𝑑𝑡. (22) Furthermore, since absorption must happen by 𝑡=∞, we know that ^∞₀𝜙(𝑝,𝑡)𝑑𝑡=1. (23) Performing integration by parts, we get the following: −1=−[𝑡𝜙(𝑝,𝑡)]^∞₀+^∞₀𝑡𝜕𝜙(𝑝,𝑡)𝜕𝑡𝑑𝑡. (24) Equations (14), (19), and (22) show that 𝜙(𝑝,𝑡) satisfies the following equation: 𝜕𝜙(𝑝,𝑡)𝜕𝜕𝑡=𝑏(𝑝)²𝜙(𝑝,𝑡)2𝜕𝑝²𝜕+𝑎(𝑝)𝜙(𝑝,𝑡).𝜕𝑝 (25) Using (25) and the fact that 𝜙(𝑝,𝑡) approaches 0 faster than t approaches ∞, we can rewrite (24) as −1=0+^∞₀𝑡𝜕𝑏(𝑝)²𝜙(𝑝,𝑡)2^∗𝜕𝑝²𝜕+𝑎(𝑝)𝜙(𝑝,𝑡)𝜕𝑝𝑑𝑡. (26) After interchanging the order of integration and differentiation we get 𝑑−1=𝑏(𝑝)²𝑡(𝑝)2^∗𝑑𝑝²𝑑+𝑎(𝑝)^̄‌𝑡(𝑝),𝑑𝑝 (27) where =𝑡(𝑝)=meantimeuntilabsorption_∞0𝑡𝜙(𝑝,𝑡)𝑑𝑡=¹₀𝑡(𝑝,𝑥)𝑑𝑥 (28) and 𝑡(𝑝,𝑥)𝑑𝑥 is the mean time that the daf spends in the interval (𝑥,𝑥+𝛿𝑥) before absorption occurs.

We are interested in the case, where 𝑝=1/2𝑁, since this is the initial frequency of the derived allele. In this case, we are interested only in values of x greater than 1/2𝑁, and for these values we can write 𝑡(𝑝,𝑥)=2𝑃₁(𝑝)∫¹_𝑥𝜓(𝑦)𝑑𝑦,𝑏(𝑥)𝜓(𝑥) (29) and 𝜓(𝑥) is defined in (17).

Under the simplest genetic model that assumes no selection and no recurrent mutation, we can set 𝑠=0 in (17) and (21) and show that 𝑃₁(𝑝) reduces to p and 𝜓(𝑦) reduces to 1. It follows from this that (29) can be reduced to 𝑡(𝑝,𝑥)=2𝑝(1−𝑥)=𝑥(1−𝑥)/2𝑁4𝑁𝑝𝑥. (30) Under a genetic model where 𝑠≠0, using 𝛾=2𝑁𝑠, (29) can be rewritten as 2𝑡(𝑝,𝑥)=2𝑁1−𝑒^−2𝛾𝑝/1−𝑒^−2𝛾∫¹_𝑥𝑒^−2𝛾𝑦𝑑𝑦𝑥(1−𝑥)(𝑒^−2𝛾𝑥. (31) After integrating and simplifying the terms, we obtain 𝑡(𝑝,𝑥)=2𝑁1−𝑒^−2𝛾𝑝1−𝑒^{−2𝛾(1−𝑥)}𝛾1−𝑒^−2𝛾.𝑥(1−𝑥) (32) Finally, substituting 𝛾=2𝑁𝑠 and 𝑝=1/2𝑁, and invoking the approximation 𝑒^−𝑎=(1−𝑎) for small values of a, 𝑡(𝑝,𝑥) reduces approximately to 2𝑓(𝑥)=𝑡(𝑝,𝑥)≈1−𝑒^{−2𝛾(1−𝑥)}1−𝑒^−2𝛾,𝑥(1−𝑥) (33) 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 ^𝑥2_𝑥1𝜃𝑓(𝑥)𝑑𝑥=^𝑥2_𝑥1𝑔(𝑥)𝑑𝑥 (34) is the expected number of sites in the population with derived allele frequency between x₁ and 𝑥₂ (where 𝜃 equals 2𝑁𝜇, 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. 2𝑔(𝑥)=𝜃1−𝑒^{−2𝛾(1−𝑥)}1−𝑒^−2𝛾𝑥(1−𝑥)=4𝑁𝜇1−𝑒^{−2𝛾(1−𝑥)}1−𝑒^−2𝛾𝑥(1−𝑥). (35) In a sample of size n, the expected number of sites with i (which ranges from 1 to 𝑛−1) copies of the derived allele is defined as a function of 𝑔(𝑥): 𝐹(𝑖)=¹₀𝑔(𝑥)𝑃(𝑖∣𝑥)𝑑𝑥=¹₀𝑛𝑖𝑔(𝑥)𝑥^𝑖(1−𝑥)^𝑛−𝑖𝑑𝑥. (36) 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 𝑋=(𝑋₁,𝑋₂,𝑋₃,…,𝑋_𝑛−1), 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 X_i is an independent Poisson distribution with mean equal to F(i) [2]. This framework allows us to define the probability of observing x_i sites that have i copies of the derived allele (and 𝑛−𝑖 copies of the ancestral allele) as the following: 𝑃𝑋_𝑖=𝑥_𝑖=𝑒∣𝜃,𝛾^{−𝐹(𝑖)}𝐹(𝑖)^𝑥𝑖𝑥_𝑖!. (37) Since the X_i’s are independent, the probability of observing 𝑋=(𝑋₁,𝑋₂,𝑋₃,…,𝑋_𝑛−1) is given as 𝑃(𝑋)=𝐿(𝜃,𝛾)=^𝑛−1_𝑖=1𝑃𝑋_𝑖=𝑥_𝑖∣𝜃,𝛾. (38) 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 2lnΛ 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.

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 (𝑠=0), the expected number of polymorphic sites with population daf x can be computed by taking the product of the per-locus mutation rate (𝜃=2𝑁𝜇) and the probability under a neutral model of a single mutation having a frequency of x (from (30)): 𝑔_neutral(𝑥)=𝜃4𝑁𝑝𝑥=2𝑁𝜇4𝑁(1/2𝑁)𝑥=4𝑁𝜇𝑥=2𝜃𝑥. (39) 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 𝑃_𝑚(𝑥)=1−(all𝑚arederived)−(all𝑚areancestral)=1−𝑥^𝑚−(1−𝑥)^𝑚. (40) 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 (𝑔_neutral(𝑥)) 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 𝑔_neutral(𝑥) and 𝑃_𝑚(𝑥) over the range of possible values for x: 𝐿(𝑚)=¹₀𝑔_neutral(𝑥)𝑃_𝑚(𝑥)𝑑𝑥=2𝜃¹₀1−𝑥^𝑚−(1−𝑥)^𝑚𝑥𝑑𝑥=2𝜃^𝑚−1_𝑘=11𝑘. (41) Finally, the total number of expected synonymous polymorphic sites in both species’ sample data is given as SP=𝐿(𝑚)+𝐿(𝑛). (42)

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 P_m(x) from 0 to 1: 𝐻(𝑚)=¹₀𝑔(𝑥)𝑃_𝑚=(𝑥)𝑑𝑥₁₀𝑔(𝑥)1−𝑥^𝑚−(1−𝑥)^𝑚𝑑𝑥. (43) Finally, the total expected number of replacement polymorphic sites in both species’ sample data is given as RP=𝐻(𝑚)+𝐻(𝑛). (44)

5.3. Expected Number of Synonymous Fixed Substitutions

When 𝑠=0, the expected number of fixed substitutions in one species relative to another that diverged 𝑡_div2𝑁 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 𝜃𝑡_div2𝑁. (45) 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 𝑒^−𝑎=(1−𝑎) for small values of a. Thus, for a newly derived neutral allele that has an initial frequency of 1/2𝑁, the probability of fixation is also 1/2𝑁.

Therefore, the total expected number of fixed substitutions in species 1 is (𝜃𝑡_div12𝑁)(2𝑁)=𝜃𝑡_div. (46) 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 ¹₀𝑇_𝑚(𝑥)𝑔_neutral(𝑥)𝑑𝑥=𝜃¹₀𝑥^𝑚2𝑥𝑥𝑑𝑥=𝜃^𝑚𝑚¹₀2=𝜃𝑚, (47) where 𝑇_𝑚(𝑥)=𝑃𝑟(aderivedalleledaf𝑥<1isobservedwith𝑥=1inasize𝑚sample) and 𝑔_neutral(𝑥) is given in (39).

Therefore, the total expected number of synonymous fixed substitutions in both species’ sample data is given as SF=𝜃𝑡_div+2𝑚+𝜃𝑡_div+2𝑛=2𝜃𝑡_div+1𝑚+1𝑛. (48)

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 𝑡_div2𝑁 generations ago is given as the product of (45) and (21): (𝜃𝑡_div2𝑁)1−𝑒^{−4𝑁𝑠𝑝}1−𝑒^−4𝑁𝑠. (49) Substituting 1/2𝑁 for p and invoking the approximation that 𝑒^−𝑎=(1−𝑎) for small values of a, we arrive at the following: (𝜃𝑡_div2𝑁)2𝑠1−𝑒^−2𝛾=𝜃𝑡_div2𝛾1−𝑒^−2𝛾. (50) 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 𝑄(𝑚)=¹₀𝑇_𝑚(𝑥)𝑔(𝑥)𝑑𝑥=2𝜃¹₀𝑥^𝑚−11−𝑒^{−2𝛾(1−𝑥)}(1−𝑒^−2𝛾)(1−𝑥)𝑑𝑥. (51) Therefore, the total expected number of replacement fixed substitutions in both species’ sample data is given as RF=𝜃2𝛾𝑡_div1−𝑒^−2𝛾+2𝐺(𝑚)+𝜃2𝛾𝑡_div1−𝑒^−2𝛾+2𝐺(𝑛)=2𝜃2𝛾𝑡_𝑑𝑖𝑣1−𝑒^−2𝛾+𝐺(𝑚)+𝐺(𝑛),where𝐺(𝑚)=𝑄(𝑚)/2𝜃. (52)

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.