Copyright Ó 2008 by the Genetics Society of America
DOI: 10.1534/genetics.108.087361

The Polymorphism Frequency Spectrum of Finitely Many Sites
Under Selection
Michael M. Desai* and Joshua B. Plotkin†,1*Lewis–Sigler Institute for Integrative Genomics, Princeton University, Princeton, New Jersey 08544 and †Department of Biology andProgram in Applied Mathematics and Computation Science, University of Pennsylvania, Philadelphia, Pennsylvania 19104
Manuscript received January 21, 2008
Accepted for publication October 6, 2008
ABSTRACTThe distribution of genetic polymorphisms in a population contains information about evolutionary
processes. The Poisson random field (PRF) model uses the polymorphism frequency spectrum to infer the
mutation rate and the strength of directional selection. The PRF model relies on an infinite-sites
approximation that is reasonable for most eukaryotic populations, but that becomes problematic when u islarge (u * 0.05). Here, we show that at large mutation rates characteristic of microbes and viruses the infinitesites approximation of the PRF model induces systematic biases that lead it to underestimate negative
selection pressures and mutation rates and erroneously infer positive selection. We introduce two new
methods that extend our ability to infer selection pressures and mutation rates at large u: a finite-site
modification of the PRF model and a new technique based on diffusion theory. Our methods can be used toinfer not only a ‘‘weighted average’’ of selection pressures acting on a gene sequence, but also the distributionof selection pressures across sites. We evaluate the accuracy of our methods, as well that of the original PRF
approach, by comparison with Wright–Fisher simulations.

M

UTATION rates and selective pressures are ofcentral importance to evolution. The number
and frequency distribution of genetic polymorphisms
within a population carry information about thesefundamental processes. Polymorphisms at higher frequencies reflect weaker selective pressures (or positive
selection), and vice versa. Similarly, a larger number ofpolymorphisms indicates a higher mutation rate. Thus
we can use the polymorphism frequency spectrum
observed in genetic sequences sampled from a population to infer the mutation rate and the selection pressure
acting on the sequence.
This intuition can be formalized into a rigorous
method for estimating selection pressures and mutation
rates by calculating the likelihood of sampled polymorphism data as a function of these parameters. The
Poisson random field (PRF) model provides an important and widely used method of doing so. The PRF
model assumes a panmictic population of constant size,
free recombination, infinite sites, no dominance or
epistasis, and equal selection pressures at all sites.
Under these assumptions, Sawyer and Hartl (1992)
showed that the distribution of frequencies of mutant
lineages in a population forms a Poisson random field
whose properties depend on the selection pressure andthe mutation rate. Hartl et al. (1994) and Bustamante

1
Corresponding author: Department of Biology, University of Pennsylvania, Philadelphia, PA 19104. E-mail: jplotkin@sas.upenn.edu

Genetics 180: 2175–2191 (December 2008)

et al. (2001) developed a maximum-likelihood methodof estimating these parameters from data on the polymorphism frequency spectrum. This method has been
widely used to study, for example, purifying selection onsynonymous (Hartl et al. 1994; Akashi and Schaeffer
1997; Akashi 1999) and nonsynonymous (Hartl et al.
1994; Akashi 1999) variation, and the evolution of basecomposition (Lercher et al. 2002; Galtier et al. 2006).
Closely related to these analyses of polymorphism
data are methods that calculate, on the basis of the PRF
model, the ratio of the expected number of polymorphisms within species to divergence between species forsynonymous and nonsynonymous sites [using the idea
behind the McDonald–Kreitman test (McDonald andKreitman 1991)]. These methods discard some of theavailable data, as they depend only on the number ofpolymorphisms and not on their full frequency spectrum. However, they are also less sensitive to assumptions (Sawyer and Hartl 1992; Loewe et al. 2006).
Such methods have been applied to estimate selection
pressures on synonymous variation (Akashi 1995), on
nonsynonymous mutations in mitochondrial genomes(Nachman 1998; Rand and Kann 1998; Weinreich and
Rand 2000), and on nonsynonymous variation in avariety of nuclear genomes (Bustamante et al. 2002;
Sawyer et al. 2003; Bartolome et al. 2005), including
humans (Bustamante et al. 2005).
Recent theoretical work has focused on relaxing
various assumptions of the original PRF method. These
include allowing for dominance (Williamson et al.

2176

M. M. Desai and J. B. Plotkin2004), population subdivision (Wakeley 2003), changing population size (Williamson et al. 2005), andlinkage between sites (Zhu and Bustamante 2005).
Several methods for studying the properties of the distribution of selection pressures across sites based on the
PRF model have also been developed, using the polymorphism frequency spectrum (Nielsen et al. 2005),
the ratio of polymorphism to divergence (Sawyer et al.
2003; Loewe et al. 2006), or several of these methods inconjunction (Bustamante et al. 2003; Piganeau andEyre-Walker 2003; Boyko et al. 2008).
A fundamental assumption of the PRF approach isthat mutation is irreversible—the infinite-sites assumption. Thus any particular site can be only transiently
polymorphic, and there is no steady-state solution to theevolutionary dynamics at any given site. However, since
the constant creation of new polymorphic sites is balanced by older polymorphisms fixing or going extinct,
the frequency distribution of polymorphisms across sites
currently polymorphic does reach an equilibrium. Thisapproach has important advantages. Since it describes
the frequency spectrum of new mutations (‘‘derived
alleles’’)—the frequency spectrum relative to the lastancestral state—it provides a natural framework foranalyzing polymorphism data when the ancestral stateis known. This ancestral state can often be inferred froman appropriate outgroup.
Despite its advantages, the infinite-sites approximation can present problems. The approximation isreasonable for most data currently available from typicaleukaryotic populations. However, in many biologicallyreasonable parameter regimes, particularly those relevant to bacterial and viral populations, more than one
mutational event may contribute to polymorphism at agiven site. In this article, we show that under theseparameter regimes the infinite-sites assumption causes
the PRF method to underestimate negative selection
pressures and mutation rates by as much as an order ofmagnitude. In addition, the PRF method often infers
that a gene is under strong positive selection when in fact
the gene is experiencing weak negative selection. This
problem arises both for inferences based on the polymorphism frequency spectrum and for inferences basedon the ratio of within-species polymorphism to betweenspecies divergence, but here we focus exclusively on theformer.
In this article, we present two methods that relax theinfinite-sites assumption of the PRF method, each with
its own advantages and drawbacks. Rather than studying
mutant lineages across a sequence, our methods focus
on explicit models of the evolutionary dynamics atindividual sites. We first present a modification of the
PRF method that retains the essential framework, but
calculates the frequency distribution of mutant lineages
at each site rather than across the whole sequence. We
next present an alternative method based on wellknown diffusion equations in place of the PRF frame-work. This alternative framework avoids all of the finitesite biases of the PRF, but it cannot make use ofknowledge about the ancestral state. Rather than describing a steady-state frequency distribution of derived
nucleotides relative to this ancestral state, it describes
the frequency distribution of all four nucleotides
possible at each site—a fundamentally different steadystate. Both of our methods allow us to estimate theselection pressure and the mutation rate from data on
the polymorphism frequency spectrum. In addition,these methods also allow us to infer the distribution of
selection pressures across sites.
To assess the accuracy of these methods, we generate
polymorphism data from simulated Wright–Fisher populations with known selection pressures and mutation
rates. By comparing inferences drawn from these simulated data sets, we demonstrate that our methods extend
and improve upon the original PRF approach. Throughout this article, we focus on accurately inferring the signand strength of negative selection, since the most troubling bias in the original PRF method is erroneous
inference of positive selection when mutation rates arelarge. We focus primarily on situations when theancestral state at each site cannot be reliably inferred,
which is the typical situation when the mutation rate islarge.
We emphasize that the effects of finite sites are ofpractical relevance only when the mutation rate isrelatively large (u per site $0.05). As a result, themethods of inference developed here are not necessary
for analyzing human or Drosophila population-genetic
data. However, as we shall demonstrate, finite-site effects
have significant practical implications when studying
the population genetics of viruses, microbes, and somehigher eukaryotes, such as sea squirts and starfish, thatexperience large mutation rates (Drake et al. 1998;
Lynch and Conery 2003).
THEORYThe Poisson random field model of polymorphisms:
We begin by outlining the PRF model of the sitefrequency spectrum developed by Sawyer and Hartl(Sawyer and Hartl 1992; Hartl et al. 1994). Thismodel assumes that mutations occur in a population ofeffective size N at a Poisson rate Nu, where u is the persequence mutation rate, and are all subject to selection ofstrength s. The fate of each mutant lineage is modeled bya diffusion approximation to the processes of selection
and drift. When a new mutant lineage enters the population, it is assumed to arise at a site that has notpreviously experienced any mutations (the infinite-sites
assumption). Each mutant lineage is assumed to be
independent of all others (the free-recombination
assumption).
Extending earlier work by Wright (1938) andMoran (1959), Sawyer and Hartl (1992) calculatedPolymorphism Frequency Spectrum

a steady-state distribution of mutant lineage frequencies. They found that the number of lineages withfrequency between x and x 1 dx is Poisson distributedwith mean f(x)dx, wheref ðxÞ ¼ ul

1  e 2gð1xÞ
1:
2g
xð1  xÞ1eð1ÞHere g [ Ns is a measure of the strength of selection on
the mutant lineages and ul [ 2Nu is twice the population
per-sequence mutation rate. The function f(x) is referred to as a Poisson random field. In other words, thenumber of mutant lineages with frequency between x1
andÐ x2 x2 is a Poisson random variable with mean
x1 f ðxÞdx. In addition, the number of mutant lineages
with frequency in [x1, x2] is independent (as a randomvariable) from the number of mutant lineages withfrequency in [y1, y2], provided these intervals do notintersect. Note that f(x) is not integrable at x ¼ 0. This
divergence occurs because the steady state arises from abalance between new mutations constantly occurringand older lineages fixing or going extinct. Thus there isno finite, steady-state expression for the number oflineages that have fixed or gone extinct.
Hartl et al. (1994) and Bustamante et al. (2001)
used Equation 1 as the basis for maximum-likelihood(ML) estimation of the mutation rate ul and selection
pressure g from polymorphism data. They imagined
sampling n individuals from a population with thissteady-state distribution of segregating mutant lineages.
They made the infinite-sites assumption that all mutant
lineages occur at different sites, consistent with theearlier assumption that each lineage is independent.
Since the number of mutant lineages at frequency x in
the population is Poisson distributed with mean f(x)dx,
the number of sampled sites containing i mutant
nucleotides (we refer to these as i-fold mutant sites) is
Poisson distributed with mean 
ð11  e 2gð1xÞ
1n i
F ðiÞ ¼ ul
x ð1  xÞni dx:
1  e 2g xð1  xÞ i0ð2ÞThis equation leads immediately to a maximum-likelihood procedure for estimating g and ul (Bustamante
et al. 2001). A set of sequences from n sampledindividuals within a population will contain somenumber, yi, of i-fold mutant sites for 0 , i , n. The setof values y1, y2, . . . , yn1 is called the site-frequency
spectrum of the observed data. The probability of aspectrum {yi}, given g and ul, isLu ðul ; gÞ ¼

n
1
Y
i¼1

e F ði;g;ul Þ½F ði; g; ul Þyi
:
yi !

ð3Þ

For an observed spectrum {yi} in a particular data set,
one maximizes this likelihood over ul and g to estimate
the mutation rate and selection pressure.

2177

The likelihood expression above assumes we know
which nucleotide is ancestral and which nucleotide is
the mutant (‘‘derived’’) at each polymorphic site. We
refer to this situation as the unfolded case. When we do
not have this information, we cannot distinguish between an i-fold mutant site and an (n  i)-fold mutant
site. In this case, a data set will contain some number, yi,i-fold and/or (n  i)-fold mutant sites, where i runs
between 1 and the largest integer #n/2. We refer to thisas the folded case. In this situation, the likelihood of aparticular data set Lf ðul ; gÞ is given by the same
expression as for Lu, but with the product running from
1 to n/2 and with F(i) replaced by F(i) 1 F(n  i) (except
if i ¼ n/2).
The infinite-sites approximation: The PRF model
makes two key assumptions: that each site is independent of all the others and that two mutant lineages neversegregate at the same site. The former assumption isequivalent to assuming free recombination between all
segregrating sites, and it has been investigated elsewhere(Akashi and Schaeffer 1997; Bustamante et al. 2001);
we return to it in the discussion. The main focus of thisarticle, however, is on the second assumption of the PRF,
that there are an infinite number of sites, which has notbeen discussed much in the literature. This assumption
can lead to problems that are most apparent when
considering how the PRF method treats ‘‘multiply polymorphic’’ sites—those that exhibit more than two typesof segregating nucleotides. Polymorphisms of this typeare indeed observed in data analyzed by the PRF method(Hartl et al. 1994). We refer to the configuration of aparticular site as (a, b, c, d), where a, b, c, and d are thenumbers of sampled sequences that exhibit each of
the four nucleotides. When we have unfolded data, a isthe frequency of the ancestral nucleotide and b, c, and d
are the frequencies of the three possible mutant nucleotides, in order of decreasing frequency. For folded
data, a, b, c, and d are the frequencies of all four possible
nucleotides, again in order of decreasing frequency. In
the original PRF analysis, a site with a (12, 1, 1, 0)
configuration, for example, is treated identically to a sitewith a (12, 2, 0, 0) configuration. Such a treatment isincorrect: the former configuration can arise only fromtwo mutant lineages, whereas the latter configuration
could be caused by a single mutant lineage (presumably
at relatively high frequency in the population). Yet the
PRF analysis excludes the first possibility and treats bothconfigurations as if they were (12, 2, 0, 0) sites (Hartl
et al. 1994; Bustamante et al. 2001). Similarly, the PRF
method treats (10, 2, 2, 0), (10, 3, 1, 0), and (10, 2, 1, 1)
sites as if they were in a (10, 4, 0, 0) configuration, etc.
The infinite-sites approximation also affects sites thatexhibit only two types of segregrating nucleotides.
When multiple lineages are segregating at the same
site, an i-fold and a k-fold mutant lineage sampled at thesame site can lead to an apparently (i 1 k)-fold mutant
site, if the two lineages happen to be mutations to the2178

M. M. Desai and J. B. Plotkinsame nucleotide. In other words, a (12, 2, 0, 0) site could
reflect two low-frequency mutant lineages or one
higher-frequency lineage, but the PRF method incorrectly assumes that only the latter is possible. This leads
to systematic biases in the estimates of selection andmutation obtained by the PRF method: by disregarding
the possibility that an apparently high-frequency mutant lineage is actually several lower-frequency mutant
lineages, the PRF method underestimates the mutation
rate u and the strength of negative selection jsj.
Since a mutant lineage will survive on averageOðln½1=jsjÞ generations before fixing or going extinct(Desai and Fisher 2007), and mutations arise at rate
Nm per site, the infinite-sites approximation will be valid
only when N m ln½1=jsj>1 (although for sufficiently
small samples the condition is weaker, since we may
never sample multiple lineages even though they aresegregating at the same site). This condition is sometimes violated in real populations, particularly of viruses
and microbes (e.g., Hartl et al. 1994).
We stress that in parameter regimes relevant to mosteukaryotes, including humans and Drosophila, finitesites biases are negligible. But in parameter regimes
relevant to bacteria and viruses, sites with multiple
segregrating mutations have a disproportionate weightin estimates of selection and mutation, and thus they
can lead to errors of an order of magnitude or more
(Figure 1). In such circumstances, the PRF method also
erroneously infers positive selection in many situations
where selection is actually negative (Figure 1).
A per-site Poisson random field model of polymorphisms: In this section, we describe a method that
extends the PRF framework to the case of finite sites andtakes full advantage of the information provided by thefrequencies of all possible configurations at a site. The
basic idea behind this modified approach is to recast the
PRF framework on a per-site basis. We describe the steadystate frequency distribution of mutant lineages at a givensite. From this, we can calculate the probability that asample of n individuals will contain any configuration ofmutants at that site. As in the original PRF method, we
retain the assumption of free recombination, so that the
DNA sequence is a collection of independent sites. Thus
our per-site analysis leads directly to ML estimation of
mutation rate and selection strength.
We begin by recasting the PRF expression for thesteady-state distribution of mutant lineages to describe
the frequencies of mutant lineages at a single site. At agiven site, we havef ðxÞ ¼ us

1  e 2gð1xÞ
1
;
2g
xð1
 xÞ
1e

ð4Þ

where us is the per-site value, us ¼ 2Nm, where m is the persite mutation rate. Using this formula to describe
multiple lineages at a single site is somewhat peculiar,
because this result assumes that all mutant lineages

behave independently of one another. Clearly this is notstrictly true, since the mutant lineages are segregating at
the same site. However, provided two mutant lineages
rarely achieve simultaneous high frequencies in the
population, then the assumption of independent mutant lineages is a good approximation. This assumptionof noninteracting mutant lineages will often hold evenwhen the other aspects of the infinite-sites approximation are violated.
Analogous to the original PRF method, at a single site
the number of mutant lineages that are observed i times
in a sample of n sequences (‘‘i-fold mutant lineages’’) is
Poisson distributed with mean 
ð11  e 2gð1xÞ1
n i
x ð1  xÞni dx:
F ðiÞ ¼ us
2g
i
xð1

xÞ
1

e
0
ð5Þ
On the basis of this, we can calculate the probability ofany particular polymorphism configuration at a site.
We begin by describing this calculation in the unfolded case. The probability that a site is monomorphic
is just the probability that no i-fold mutant lineages arefound at that site, for all i between 1 and n  1. This is
Pmono ¼ Pðn;0;0;0Þ ¼ e F ð1Þ e F ð2Þ . . . e F ðn1Þ :ð6Þ

We define M [ Pmono to denote this monomorphism
probability. The probability that a site exhibits a (n  1,
1, 0, 0) configuration is the probability that a single 1fold mutant lineage is sampled, and no 2-fold or higher
lineages are found,Pðn1;1;0;0Þ ¼ F ð1ÞM ;ð7Þwhere M denotes the probability of monomorphism as
above. The probability that a site exhibits an (n  2, 2, 0,0) configuration is more complex. This configuration
could arise from a single 2-fold sampled lineage (asassumed under the standard PRF method) or it could
arise from two 1-fold sampled mutant lineages that
happen to involve mutations to the same nucleotide.
Hence its probability isPðn2;2;0;0Þ ¼

F ð1Þ2 1M 1 F ð2ÞM ;
3
2!

ð8Þ

where the factor of 13 is the probability that two
mutations result in the same nucleotide. This expression assumes that mutations between all possible nucleotides are equally likely—the obvious generalization
applies when there are mutational biases, which we do
not discuss further. Similarly the probability of an (n  2,
1, 1, 0) configuration isPðn2;1;1;0Þ ¼

F ð1Þ2 2
M :
3
2!

ð9Þ

The probabilities of more complex configurations
can be calculated in a similar way. The probability ofPolymorphism Frequency Spectrum

an (n  4, 4, 0, 0) configuration, for example, is theprobability of four 1-fold lineages to the same nucleotide plus the probability of two 1-fold lineages and a 2fold lineage, plus the probability of two 2-fold lineages,
plus the probability of one 1-fold lineage and one 3-foldlineage, plus the probability of a single 4-fold lineage.
We haveF ð1Þ4 1F ð1Þ2 F ð2Þ 1
F ð2Þ2 1Pðn4;4;0;0Þ ¼ M
1
1
2!
4! 33
2! 3
32

1
ð10Þ
1 F ð1ÞF ð3Þ 1 F ð4Þ :3
In general, the probability of a particular configuration
is given by the sum of the probabilities of all possible
partitions of n that lead to that configuration.
In the folded case, these calculations become fairly
complex. The probability of a (12, 2, 0, 0) folded
configuration, for example, includes the probability ofa single 12-fold sampled lineage, as well as two 6-foldsampled lineages, and so on. Thousands of terms may
arise in the expression for the probability of a particularconfiguration, even for moderate values of n. We do not
quote any of these expressions here, but rather we havedeveloped a computer program to output symbolic expressions for the probabilities of all possible folded aswell as unfolded configurations, for a given sample sizen. The algorithm used in this program is described in the
appendix, and the program is freely available on
request.
These probabilities of site configurations form thebasis of maximum-likelihood parameter estimation.
The probability of a data set with L total sites, including
ya,b,c,d sites in configuration (a, b, c, d), is given byY yL!
a;b;c;dQ
Pa;b;c;d;ð11Þ
ya;b;c;d !
where the products are over all possible configurations.
Given a particular data set, we maximize this probability
over us and g to find the ML estimate of theseparameters. In the original PRF method, this maximization is particularly simple, because the ML estimate
for u can be expressed analytically in terms of the ML
estimate for g, leaving a one-dimensional numericalmaximization procedure to estimate g. In our per-site
PRF method, however, a full two-dimensional numericalmaximization is required to find the ML estimates of us
and g. We have implemented a computer program to
perform this maximization; it is freely available on
request.
Our per-site PRF method relaxes most of the consequences of the infinite-sites approximation inherent
in the original PRF estimation procedure. We allow for
the possibility that multiple mutant lineages contributeto the polymorphism observed at a single site. Thus weavoid the systematic underestimation of u and g, andincorrect inference of positive selection, that affect the2179

traditional PRF when mutation rates are large (Figure1). This new method also uses all of the data available inthe sample, including the differences between (n  2, 2,0, 0) and (n  2, 1, 1, 0) sites and the number ofmonomorphic sites (note the infinite-sites approximation means that the original PRF cannot predict the
number of monomorphic sites and hence makes no useof these data). It does still retain one aspect of theinfinite-sites approximation: it assumes that mutant
lineages segregating at the same site do not interferewith each other. This is never strictly true, but it is a good
approximation unless multiple mutant lineages reachhigh frequency at a given site at the same time. Note that
because of this assumption, the probabilities of all
possible configurations (a, b, c, d) described above do
not precisely sum to unity, because our approach allows
a typically small but nonzero probability of multiple
mutant processes adding to more than n sampledindividuals. The no-interference assumption is substantially weaker than the infinite-sites assumption, and thus
our revised sampling method extends the applicabilityof the PRF framework.
A per-site diffusion model of polymorphisms: In thissection, we describe a method that shifts fundamentallyfrom the PRF framework and avoids all of the problems
associated with the infinite-sites assumption. Rather
than studying the distribution of the frequencies of mutant lineages, we focus on the evolutionary dynamics ateach individual site, without keeping track of individual
mutant lineages. We develop this into a maximum-likelihood estimation of g and u from polymorphism data,
which requires neither the infinite-sites nor the nointerference approximation described above. As in theoriginal PRF method, we assume free recombination
between sites.
At an individual site, we imagine that one nucleotide ispreferred and the other three have the same fitness
disadvantage s (s , 0). We assume that mutations occur
at rate m and hence at rate m/3 between any two specificnucleotides (i.e., no mutational biases). These assumptions simplify the discussion, but are not essential. In
fact, one advantage of this approach is that these assumptions can be easily relaxed with obvious generalizations (noted below).
We can analyze the process of mutation, selection, anddrift at a single site with a three-dimensional diffusion
approximation and calculate the joint steady-state probability distribution of the frequencies of the four
possible nucleotides at the site. This then leads naturally
to the likelihood of any configuration of polymorphism
data at the site as a function of g and us and hence to ML
estimation of these parameters from data. Before exploring this fully, we first analyze a simplification of thismodel, in which we treat all three disfavored nucleotides
as a single class. Such a treatment reduces to a standard
one-dimensional diffusion process whose steady-state
probability distribution describes the frequency of the2180

M. M. Desai and J. B. Plotkin

preferred nucleotide vs. the sum of the frequencies ofthe disfavored ones; this treatment is essentially a steadystate version of Williamson et al. (2005). This simplified
model is not likely to be particularly useful, because it
requires us to know a priori the identity of the preferred
nucleotide, which in practice is part of what we wish toinfer. This simplification also discards some of theinformation in the data [e.g., not making use of thedifference between (12, 1, 1, 0) and (12, 2, 0, 0) sites].
However, first studying this analytically and computationally simpler model makes the analysis of the fullmodel more clear. We therefore begin by describing the
one-dimensional method and then turn to the threedimensional method.
The one-dimensional diffusion model: We begin bydescribing a simplified diffusion approach that calculates the frequency distribution of favored vs. disfavorednucleotides, similar in spirit to that introduced by
Mustonen and Lassig (2007). As noted above, we forsimplicity assume that one nucleotide is preferred andthe other three nucleotides are disfavored. We denotethe sum of the frequencies of the three disfavoredalleles by x; the frequency of the preferred nucleotide is
1  x.
We assume that mutation, selection, and random drift
occur at each site according to standard Wright–Fisher
dynamics. Thus the probability distribution of x can be
described by the diffusion equation
2

@
1 @
@
½mðxÞf ðx; tÞ; ð12Þ
f ðx; tÞ ¼
2 ½vðxÞf ðx; tÞ 
@t
2 @x
@x
where f(x, t) is the probability that the disfavorednucleotides have frequency x at time t and m(x) andv(x) are given bymmðxÞ ¼ sxð1  xÞ 1 mð1  xÞ  x
ð13Þ
3vðxÞ ¼

xð1  xÞ;N

ð14Þ

where s is the selection coefficient against the disfavorednucleotides (s , 0) and m is the per-site mutation rate
per individual per generation. This diffusion equation iswell known (Ewens 2004) and has the steady-statesolution derived by a zero-flux boundary condition:f ðxÞ ¼ Cx us 1 ð1  xÞus =31 e 2gx :ð15ÞHere C is a (us- and g-dependent) normalization factor,
and as before us [ 2Nm and g [ Ns.
If the frequency of disfavored nucleotides at a site
equals x, the probability that we find i such nucleotides
in a sample of n individuals is ð ni Þx i ð1  xÞni . Averaging
over x, the overall probability that we sample i disfavorednucleotides at a given site is  ð1n
F ðiÞ ¼
Cx us 1i1 ð1  xÞus =31ni1 e 2gx dx: ð16Þi0This integral, including the calculation of the normalization factor C, can be solved analytically. We find
F ðiÞ
¼

 
n Gðni 1 us =3ÞGði 1 us Þ1 F1 ði 1 us ; n 1 4us =3; 2gÞ;Gðus =3ÞGðus Þ1 F1 ðus ; 4us =3; 2gÞið17Þ
where G is Euler’s Gamma function and 1F1 is a hypergeometric function.
The expression above leads immediately to a maximum-likelihood method for estimating g and us in theunfolded case. Note, however, that in contrast to theoriginal PRF method, this model describes a steady stateof the frequency distribution of preferred relative tounpreferred nucleotides, rather than derived vs. ancestral nucleotides. This is a very different sort of steadystate, a point to which we return in more detail below.
Since this model refers only to the preferred andunpreferred states at each site, it makes sense to apply
to unfolded data only if we assume that the ancestral
state is preferred. In a sample of n sequences each oflength L, we count the number of sites at which i
disfavored nucleotides are sampled, yi, for 0 # i # n.
Since all sites are assumed independent, each with thepolymorphism frequency distribution described above,the likelihood of the data given the parameters is
L!
Lu ðus ; gÞ ¼ Qn

nYi¼0 yi ! i¼0

F ðiÞyi :ð18Þ

For any set of polymorphism data, it is straightforward to
maximize this function numerically, producing ML
estimates of us and g. As with the per-site PRF method,
this procedure involves a two-dimensional maximization routine. Simulated data and fits confirm that thismethod produces accurate and unbiased estimates ofmutation rate and selection strength (data not shown).
When we do not know which nucleotide is preferred
at each site, which will be typical in most microbial
applications, we must use the ‘‘folded’’ version of thedata. This presents difficulties. Imagine a site with an (a,b, c, d) polymorphism configuration. We might naively
suppose that since any of the four nucleotides could be
the preferred one, the probability of these data is simply
F(b 1 c 1 d) 1 F(a 1 c 1 d) 1 F(a 1 b 1 d) 1 F(a 1 b 1 c).
However, this is not the case. For example, if a is indeed
the preferred nucleotide, then F(b 1 c 1 d) does notequal the probability that the three disfavored nucleotides will form a (b, c, d) configuration. Rather, it equals
the sum of the probabilities that the three disfavorednucleotides will form a configuration (i, j, k), summed
over all i, j, k triplets that sum to b 1 c 1 d. This is a serious
problem, because the difference between F(b 1 c 1 d)and the probability of the data depends on b 1 c 1 d,and hence it is not identical for all four possible preferred nucleotides. Simply assuming that the most common

Polymorphism Frequency Spectrum

nucleotide is the preferred one (Hartl et al. 1994) is areasonable approach. But this approach will be inaccurate for sites where the most common nucleotide
is not overwhelmingly so; when this situation describes a
substantial fraction of sites, the method will fail. As aresult, the one-dimensional diffusion framework does
not allow for a rigorous ML estimate of parameters withfolded data.
To perform rigorous ML fits to folded frequency data
we must turn to a three-dimensional diffusion method.
The three-dimensional diffusion model: Rather thanconsidering all disfavored nucleotides as a single class,
we can instead keep track of the evolutionary dynamics
of all four possible nucleotides at a site. To do so, we
assume the standard four-allele Wright–Fisher dynamics, with mutation at rate m=3 between any two particular
nucleotides, and selection acting with strength s against
the three disfavored nucleotides. The dynamics canthen be described by a three-dimensional diffusion
equation for the joint distribution of the frequenciesof the three disfavored alleles x1, x2, and x3, f(x1, x2, x3, t)(where the preferred allele has frequency x0 ¼ 1  x1 
x2  x3). We have@f ðx1 ; x2 ; x3 ; tÞ
@t
3 X
3
3
 X
1X
@2 
@
½ Mi 3 f ;
Vij 3 f ¼
2 i¼1 j¼1 @xi @xj
@xi
i¼1
where
Mi ¼ ð1 1 sÞxi 1 

3
X
j¼1

Vii ¼

!
xj

1

mð1  4xi Þ3xi ð1  xi ÞN

Vij ¼ 

xi xj
;
N

ð19Þ

ð20Þ

ð21Þ

2181

Averaging over f, we therefore find that the probabilityof sampling a site in an unordered (n0, n1, n2, n3)configuration isn!
Pn0 ;n1 ;n2 ;n3 ¼ C
n0 !n1 !n2 !n3 !

ð 1 ð 1x1 ð 1x1 x2
0

0

ð25Þ

0

3 e 2gðx1 1x2 1x3 Þ x1z1n1 x2z1n2 x3z1n3 ½1  x1  x2  x3 z1n0
3 dx3 dx2 dx1 :ð26Þ
If we have unfolded data and wish to assume that theancestral state is preferred, we can use unfolded ML
inference. In this case, the probability of finding a site inan ordered unfolded configuration (a, b, c, d) (where by
convention a is the number of individuals that have thepreferred nucleotide and b $ c $ d) isX
uPa;b;c;d¼Pn0 ;n1 ;n2 ;n3 ;ð27Þ
fn0 ;n1 ;n2 ;n3 gwhere the sum is over all unordered configurations (n0,n1, n2, n3) that give rise to the ordered unfolded
configuration (a, b, c, d).
If on the other hand we do not have a reliable
outgroup (or choose to ignore the information froman outgroup on the ancestral state), we can use folded
ML inference. Here, the probability of sampling a site in
the ordered configuration (a, b, c, d) is just the sum of
the probabilities assuming that each of the four possible
nucleotides is preferred. As before, we adopt theconvention that ordered folded configurations are
written as (a, b, c, d) with a $ b $ c $ d. The probabilityfof a folded configuration Pa;b;c;dis thenXf¼
Pn0 ;n1 ;n2 ;n3 ;ð28Þ
Pa;b;c;dfn0 ;n1 ;n2 ;n3 gði 6¼ jÞ:

ð22ÞThis is a somewhat less well-known diffusion equation(Wright 1949; Watterson 1977); the steady-statesolution isf ðx1 ; x2 ; x3 Þ ¼ C½x1 x2 x3 ð1  x1  x2  x3 Þz e 2gðx1 1x2 1x3 Þ ;ð23Þ
where C is a normalization factor, and we have definedz ¼ us/3  1.
Given the frequencies x0, x1, x2, and x3 of nucleotides
in the population, the probability of sampling a site in
an unordered configuration (n0, n1, n2, n3) in a sample ofn individuals (adopting the convention that the firstnucleotide listed is the preferred one) is just themultinomial probabilityn!
ð1  x1  x2  x3 Þn0 x1n1 x2n2 x3n3 :n0 !n1 !n2 !n3 !

ð24Þ

where in this case the sum is over all unordered
configurations that give rise to the ordered folded
configuration (a, b, c, d).
For either folded or unfolded data, given n samples
of a sequence L sites long, with ya,b,c,d sites in an (a, b, c, d)polymorphism configuration, the likelihood of thedata isYL!
½Pa;b;c;d ya;b;c;d ;Lðus ; gÞ ¼ Q
ð29Þ
ya;b;c;d !
where the products are taken over all possible configurations (a, b, c, d) and Pa,b,c,d is the folded or unfolded
probability defined above. We can numerically maximize this function to find the ML estimates of us and g.
We have implemented a computer program to perform
this maximization; it is freely available on request.
Variable selection pressures across sites: Both theoriginal PRF method and the per-site methods we haveproposed in this article assume that all sites experience
the same selective pressure. In reality, we expect that

2182

M. M. Desai and J. B. Plotkin

there is some distribution of selective pressures acrosssites. Hartl et al. (1994) suggest that in this case the ML
estimate of g from their PRF method reflects a weighted
average selection pressure across the sites, but thenature of this weighting is not well understood. Almost
no weight is given to sites at which jgj?1, because thesesites will likely be monomorphic and hence ignored by
the original PRF method. It is unclear how sites withdifferent values of g of order 1 will be weighted or how
the presence of some effectively neutral sites (jgj>1)
will change the ML estimate. Piganeau and EyreWalker (2003), Nielsen et al. (2005), and (Boykoet al. 2008) have addressed these questions by introducing procedures that allow for inference of some aspectsof the distribution of selective coefficients across sites
within the PRF framework. In this section we describe asimilar generalization in our models.
The issue of variable g across sites is of particularconcern for the methods we have proposed, because
these methods make use of the monomorphism data. Ifsome number Ll of sites are effectively lethal (i.e., havejgj?1), these sites will all be monomorphic, which willtend to depress our ML estimate of us and increase ourestimate of jgj. Fortunately, our methods are able to use
the monomorphism data to investigate the number oflethal sites, Ll, or more generally the full distribution of
selection pressures across sites.
Since all of the methods we have proposed aredefined at a per-site level, it is straightforward to assume
that there are multiple different classes of sites withdifferent values of g. We can posit that there are k classesof sites. Each class is represented by Lj of L total sites andhas its own value of g, which we call gj. The probability
that a site is in an (a, b, c, d) configuration is then
Pa;b;c;d ¼
j

k
X
Lj

L
j¼1

j

Pa;b;c;d ðgj ; us Þ;

ð30Þ

where Pa;b;c;d (gj, us) is the probability that a site withparameters gj and us is in the configuration (a, b, c, d).
This expression is correct for both the per-site PRF andthe diffusion approaches.
Given our new definition of Pa,b,c,d, we can construct
the folded or unfolded likelihood of the overall polymorphism data set in exactly the same way as before.
This likelihood function now depends on 2k 1 1parameters: us, the gj, and the Lj. We can find ML
estimates of all of these parameters, using a multidimensional numerical maximization of the likelihood
function. By choosing k, we determine the resolution atwhich we measure the distribution of values of g acrosssites. Naturally, the larger the k we choose, and hence thegreater the resolution on g, the more data we require toobtain accurate estimates of the individual Lj and gj.
Rather than estimating both the Lj and the gj, we
could instead posit that there are several classes ofmutations with different prespecified gj and estimate only

the values of Lj. In other words, we ask what fraction ofsites have different values of selective constraints. We
describe here one particularly important example of a
hybrid between these two procedures, with two classes ofsites (k ¼ 2). Rather than fitting an ML estimate of g toboth classes, we assume that one class of sites is unable toevolve: mutations at these sites are lethal (more precisely, they have jgj?1). We wish to calculate thenumber of lethal sites and the average selective pressure
on the remaining, nonlethal sites. Thus we have threeparameters: the mutation rate us, the number of lethalsites L2, and the strength of selection g acting at theother L1 ¼ L  L2 sites. The probability that a site is
monomorphic is given by
Pmono ¼

L2 L1 1
1 Pmono :
L
L

ð31Þ

1
is the probability that a site with strength ofHere Pmono
selection g and mutation rate us will be monomorphic,
as defined by either the per-site PRF or the per-site
diffusion approach (whichever method we are using).
The probability that a site is in a nonmonomorphic (a, b,c, d) configuration is

Pa;b;c;d ¼

L1 1
P
:
L a;b;c;dð32Þ

Now we can write the likelihood of the data in the usual
way. This results in a three-dimensional ML problem.
However, we can simplify the problem by first maximizing L2 given g and u. We find that the ML estimate of L2
is
L̂2 ¼

1
Lmono  LPmono;
1
1  Pmono

ð33Þ

where Lmono is the number of monomorphic sites in the1data and Pmono
is the probability a nonlethal (i.e., an L1)
site is monomorphic. Substituting this value for L2, weare left with a two-dimensional maximization problem
in g and us, similar to the original situation.
It is worth exploring how this procedure for estimating the number of lethal sites utilizes the data. It turns
out that this procedure is equivalent to ignoring themonomorphism data when finding the maximum-likelihood estimates of g and us for the nonlethal sites. The
likelihood of the data ignoring monomorphic sites is"#ya;b;c;d1
Y Pa;b;c;dyp !
Lðus ; gÞ ¼ Q;
ð34Þ
1ya;b;c;d1  Pmono
where yp is the total number of nonmonomorphic sites
and the products are over all configurations of nonmonomorphic sites. After finding ML estimates of g andu from this monomorphism-ignoring likelihood function, the procedure then calculates the number ofmonomorphic sites that would be expected given g
1
1
and us. This is L1 Pmono¼ ðL  L2 ÞPmono. We thenPolymorphism Frequency Spectrum

estimate the number of lethal sites L2 as the difference
between the observed number of monomorphic sites
Lmono and the number that would be predicted if all1
sites were of the L1 variety, L̂2 ¼ Lmono  ðL  L̂2 ÞPmono.
Rearranging this expression, we see that it is identical toEquation 33 above. And indeed, plugging Equations
31–33 into Equation 30 yields Equation 34.
Thus, this procedure ignores the monomorphism
data when calculating g and us (at the nonlethal sites),
and it instead uses the monomorphism data to infer one
aspect of the distribution of g across sites—specifically,
the number of lethal sites. Since the original PRF
method also ignores monomorphism data, we obtain
this information on the distribution of g for ‘‘free,’’relative to the power of the original method, simply byshifting to the per-site model. If desired, we can also
posit that there are a number of sites L3 that are
effectively neutral (i.e., with jgj>1) and estimate L3.
This would devote some part of the data describing
polymorphisms at intermediate frequencies to estimating L3. From this procedure we could estimate thenumber of effectively lethal sites, the number ofeffectively neutral ones, and the ‘‘weighted average’’selection pressure acting on the remaining sites. If more
resolution is desired, and enough data are available, we
can increase the number of classes of sites and obtain
ML estimates of the numbers of sites in each class andthe selection pressure acting on each class.
METHODS

Simulations and fits: We used Wright–Fisher simulations to test the inferential accuracy of the PRF
method as well as the accuracy of our two alternative
methods. The Wright–Fisher model (or, more precisely,
its diffusion limit) forms the basis of the PRF method,and it is therefore the appropriate simulation framework for testing the method.
All simulations assumed a constant population of N ¼
1000 haploid individuals. Each of L ¼ 1000 sites, simulated independently, could assume one of four states: a,c, t, or g. One state is assigned fitness 1, and the otherthree states fitness 1 1 s (where s , 0). Mutations occurred at rate m per site. The allele frequencies evolved
according to the standard Wright–Fisher Markov chain(Ewens 2004). Each simulation was run for at least 10/m
generations, to ensure relaxation to steady state. At theend of the simulation, n ¼ 14 individuals were sampledfrom the population and the polymorphism frequency
spectrum was recorded. We chose to consider samples ofsize n ¼ 14 to facilitate comparison with Hartl et al.
(1994). This choice does affect the relative accuracy of
the techniques; the finite-sites biases in the PRF method
increase with n (see results). We have focused our fits
on the case of folded frequency spectra, which will
typically be the only reliable type of polymorphism data
available when us is large.

2183

We performed simulations over a wide range ofparameter values relevant to viruses and microbes. We
considered five different values of us: 0.05, 0.1, 0.5, 1.0,
5.0. For each value of us, we performed one simulation ateach of 17 different values of g, ranging from g ¼ 10.0
to g ¼ 0.1. For each set of simulation parameters (g,us), once the simulated folded polymorphism data had
been generated, ML parameter estimates ðĝ; ûs Þ wereobtained by numerical maximization of the likelihood
function, as specified by the original PRF model, theper-site PRF model, or the three-dimensional diffusion
model. This maximization can be difficult in the diffusion case (see below). A 95% confidence interval for g
was constructed according to Bustamante et al. (2001):
the interval includes those values of g within 0:5x21;0:95
log-likelihood units from ĝ (note these confidence intervals rely on the assumption of no linkage; seediscussion). The estimated parameters shown in Figure 1 are somewhat ‘‘jagged,’’ because the inference
methods have been applied to a single draw of n ¼ 14
sequences for each set of simulation parameters, as
opposed to averaging over many such draws. Figure 1
also shows fits obtained from inference techniquesapplied to data generated under the infinite-sites modelat us ¼ 0.01 (i.e., deviates drawn from a Poisson
distribution with mean given by Equation 1).
As discussed above, the original PRF model disallows
multiple mutant lineages at a site. Therefore, when a site
sampled from the simulated data exhibited more than
two types of segregating nucleotides, the frequencies of
all unpreferred nucleotides were summed to represent
the frequency of the ‘‘mutant’’ type for the purposes offitting using the original PRF model, as suggested by
Hartl et al. (1994). When fitting folded data using the
PRF method, the most common nucleotide was assumed to be the ancestral type, as suggested by Hartl
et al. (1994). This approach is not entirely accurate, as
discussed above, but it is probably the best option
available within the original PRF framework.
Numerical maximization procedure: In practice, theML estimation procedure for the diffusion model isdifficult to implement, because of the triple integral in
the definition of Pn0 ;n1 ;n2 ;n3 (as well as that implicit in thedefinition of the normalization constant C). This integral cannot be solved exactly, and it is difficult toevaluate numerically because the integrand may diverge(though the integral itself converges) near the boundary of the simplex over which it is integrated. We
adopt a hybrid method to simplify the evaluation of this
triple integral. Near the boundary of the simplex, weTaylor expand the integrand and integrate it analytically. Away from the boundary, the integrand is well
behaved and standard numerical integration has no
difficulties.
This approach is most easily achieved by making thesubstitutions y ¼ x3/(1  x1  x2) and z ¼ x2/(1  x1).
On doing so, we can rewrite the integral as2184

M. M. Desai and J. B. Plotkinð1 ð1 ð1n!
n0 !n1 !n2 !n3 ! 0 0 0
3 exp½2gðx 1 y 1 z  xy  zy  xz 1 xyzÞ

Pn0 ;n1 ;n2 ;n3 ¼ C

3 x n1 1z z n2 1z yn3 1z
3 ð1xÞn0 1n2 1n3 13z12 ð1zÞn0 1n3 12z11 ð1yÞn0 1z :

ð35ÞThis expression is much easier to handle than our
original expression, because the three integrals can bedone in arbitrary order. We divide each of the three
integrals into three pieces: one from 0 to d, one from d
to 1  d, and one from 1  d to 1. Thus the triple integral
is split into 27 total terms. For each of the integrals from0 to d or 1  d to 1, we Taylor expand the integrand in
the integration variable, and we solve the integral
analytically. All of the remaining integrals, from d to1  d, are done numerically. We must choose d large
enough that we can perform the numerical integrals
quickly, but not so large that the Taylor expansions used
for the analytical parts become invalid. For these Taylor
expansions, we need d>1=2g, d>1, d>1=ð3z 1 2Þ,
d>1=ð2z 1 1Þ, and d>1=z. For the computational
analysis described in this article, we choose whicheverof these conditions is most restrictive and set d to be onetenth of the most restrictive requirement. We find that
this choice of d is sufficiently small to provide accuracy
in the analytical parts of the integrals, but large enoughto enable quick numerical integration on the interior of
the simplex.
We have written a computer program implementing
this numerical integration and the resulting ML estimation for the diffusion model. This program and thesimpler programs for ML estimation using the onedimensional diffusion model and the per-site PRF are
freely available on request.

RESULTS

Using data from the Wright–Fisher simulations described above, we tested the inferential accuracy of theoriginal PRF method, our per-site PRF method, and thethree-dimensional diffusion method. Figure 1 compares the accuracy of selection pressures estimated fromfolded data using each of these three methods. For large
us, the original PRF method systematically underestimates the strength of negative selection, by as much as afactor of 10. In addition, the PRF method often
erroneously infers strong positive selection when in fact
mutants are under negative selection. These problems
are more severe when the mutation rate is large. The
smallest mutation rate at which such problems occur(us ¼ 0.05) is four times smaller than the mutation rateestimated for bacterial genes (Hartl et al. 1994). The
per-site version of the PRF method that we have developed improves upon the standard PRF method, but it

too exhibits systematic biases, especially when selection
is weak and the mutation rate large (Figure 1). The onedimensional diffusion method provides accurate andunbiased estimates of g over the full range of selective
pressures and mutation rates (not shown). Like its onedimensional counterpart, the three-dimensional diffusion method also provides accurate and unbiased
estimates over the full range of simulated parameters(Figure 1). When selection is weak (i.e., jgj , 1),
however, the confidence intervals on diffusion-based
estimates of g are appreciably larger. This behavior
makes perfect sense: when selection is nearly neutral
and the ancestral state is unknown, the frequency
distribution does not exhibit sufficient skew to deduce
the preferred nucleotide. As a result, the diffusionbased estimator cannot distinguish between weak positive and weak negative selection in the absence of
information on the preferred nucleotide. Thus, theconfidence intervals obtained under the folded diffusion technique properly reflect our inability to estimate
the selection pressure precisely when selection is weak.
As shown in Figure 1, when selection is weakly
negative, the original PRF method erroneously infers
positive selection. This problem occurs in the parameter
regimes that have been estimated from biological datafrom bacterial and viral populations. For example, on
the basis of n ¼ 14 sampled sequences each 367 sites
long, Hartl et al. (1994) estimated g ¼ 1.34 and us ¼
0.183 for silent sites in a bacterial gene. If we simulate
367 Wright–Fisher sites under these parameters andsample n ¼ 14 sequences, we find that the most likely
parameters fitted using the original PRF method areðĝ; ûs Þ ¼ ð118:45; 0:067Þ. That is, in this example using
estimated microbial parameter values, the PRF method
is strongly biased.
Figure 2 shows the accuracy of estimated mutation
rates using the original PRF method, the per-site PRF
method, and the three-dimensional diffusion method.
The original PRF method systematically underestimatesthe mutation rate. Estimates obtained using the per-site
PRF method are an improvement, but still exhibit somebiases. The diffusion-based method provides accurate
and unbiased estimates of the mutation rate, across thefull range of mutation rates and selection pressures(Figure 2).
In addition to the simulations described above, we
also fitted data that had been generated under theinfinite-sites model at a small mutation rate, us ¼ 0.01
(Figure 1). In this case, we generated a sampled polymorphism frequency spectrum using Poisson deviateswith mean given by Equation 1, for n ¼ 14 samples of
L ¼ 10,000 sites. When applied to such data, inferences
based on the original PRF method must be unbiased.
Figure 1 shows that inferences based on the modified
per-site PRF method and the diffusion method are also
accurate and unbiased when applied to such data. Thus,
the new methods we have developed here perform asPolymorphism Frequency Spectrum

2185

Figure 1.—Maximum-likelihood estimates of selection pressures obtained under the PRF method, the modified (per-site) PRF
method, and the three-dimensional diffusion method. Selection pressures were estimated from the folded polymorphism frequencies among n ¼ 14 sequences sampled from a simulated Wright–Fisher population. In each panel, the simulated selection
pressure g is shown on the x-axis and the estimated selection pressure ĝ on the y-axis. Dashed lines indicate 95% confidence
intervals around the ML estimates obtained from the x2-distribution (see methods). The line ĝ ¼ g is shown in red, and the lineĝ ¼ 0 is shown in gray. The diffusion method and, to a lesser extent, the modified PRF method correct the biases inherent in theoriginal PRF method. When selection is weakly negative, the diffusion method cannot reject positive selection on the basis offolded data, as indicated by the lack of the upper confidence interval for some of the fits. The simulated data in the top panels,
u ¼ 0.01, were generated according to the infinite-sites model (see methods).

well as the original PRF method under the assumptionof infinite sites and at small mutation rates.
The methods we have developed in this article also
allow us to estimate the distribution of selection pressures across sites. In one simple case discussed above, wehave presented a procedure for estimating the numberof lethal sites and the selective pressure operating on theremaining, nonlethal sites in a gene. This procedure
involves estimating g and us on the basis of polymorphic
sites alone and thereafter estimating the proportion ofobserved monomorphic sites that are lethal. To assess

the power and accuracy of this approach, Table 1 showsthe error in the predicted number of monomorphic
sites in each of our simulations, compared to thenumber of monomorphic sites actually observed. Acrossa large range of selective pressures and mutation rates,
this approach typically estimates the number of (nonlethal) monomorphic sites within a few percent. As aresult, for a gene of length L ¼ 2000 sites, at one-half ofwhich mutations are strongly deleterious (jsj?1=N ),our procedure will accurately predict the number oflethal sites within a few percent, and it will accurately

2186

M. M. Desai and J. B. Plotkin
TABLE 1
Accuracy of the estimated number of monomorphic sites
Actual us
0.05
0.1
0.5
1.0
5.0

Figure 2.—Maximum-likelihood estimates of the mutation
rate, us ¼ 2Nm, obtained under the PRF method (blue), themodified (per-site) PRF method (green), and the diffusion
methods (red). Mutation rates were estimated from thefolded polymorphism frequencies among n ¼ 14 sequences
sampled from a simulated Wright–Fisher population. The
simulated mutation rate us is shown on the x-axis and the estimated mutation rate ûs on the y-axis. The line ûs ¼ us isshown in black. For each value of us, simulations and fits
are shown for 17 different values of g, ranging from g ¼
10.0 to g ¼ 0.1. The PRF method systematically underestimates the mutation rate, especially when selection is weak.
The diffusion method provides accurate and unbiased estimates of the mutation rate across the full range of parameters.

predict the selection pressure on the remaining, nonstrongly deleterious sites.
All of the simulation results that we have presented
thus far have been for the case n ¼ 14. This choice affects
the relative accuracy of the methods. As n increases, theprobability of seeing a high-frequency polymorphism
within the sample depends more strongly on theprobability that a site has a high-frequency polymorphism in the population, f(x) for large x. This is precisely
the part of frequency distribution that is artificiallysuppressed by the infinite-sites assumption of the PRF
method. Thus, as the sample size increases the biases in
the PRF method should become more severe. We havetested this expectation by sampling n ¼ 100 simulated
Wright–Fisher sequences. In the case of us ¼ 0.1, forexample, when fitting unfolded polymorphism frequency spectra measured relative to the (known)
ancestral state at each site, the PRF method always
rejects the true value of g in favor of a weaker or evenpositive selection coefficient, throughout the range g 2[10, 0.5]. Inferential biases associated with large
sample sizes and mutation rates will become increasingly important as the availability of sequence data
improves—in metagenomic surveys of microbial populations, for example.

DISCUSSION

The Poisson random field model (Sawyer andHartl 1992) and the associated likelihood procedure

Average error (%)

Median error (%)

7.6
3.3
2.3
1.4
0.1

7.8
1.8
1.1
0.5
0.1

Accuracy of estimates for the expected number of monomorphic sites. We calculated maximum-likelihood estimatesof g and us using the one-dimensional diffusion method applied to unfolded simulated data excluding monomorphic
sites. From these values, we calculated the expected numberof monomorphic sites. Shown are the differences between theobserved and the expected number of monomorphic sites (ina simulated gene of length L ¼ 1000 sites). For each value ofus, we show both the average and the median differences in
numbers of monomorphic sites, across simulations with granging from 0.1 to 10. These results imply that we canaccurately estimate the number of sites that are monomorphic due to drift. Thus, if a gene contains a similar or larger
number of lethal sites, we can also estimate the number ofsuch lethal sites to within the above accuracy.

for estimating parameters (Hartl et al. 1994) are
perfectly valid when the assumptions underlying themethod are met—namely, infinite sites, free recombination, and constant selective pressure across sites.
Moreover, the assumption of infinite sites is perfectly
reasonable for most eukaryotic populations. It will holdwhenever N m ln½1=jsj>1; as we have seen, a rough ruleof thumb is that this tends to be true when us , 0.05. We
have shown, however, that when the mutation rate isrelatively large the PRF method can lead to incorrect
inferences. In practice, such concerns arise whenstudying populations of viruses, microbes, and thoseeukaryotes that experience large mutation rates (us .
0.05). We have developed two new methods that relax orremove the infinite-sites assumption. These new methods also extend the types of inferences that can be
drawn from polymorphism data to include inferenceson the distribution of selective pressures across sites.
It may seem surprising that the infinite-sites approximation can lead to substantial errors in the mutation
rates and selective strengths inferred by the PRF
method. After all, sites at which multiple mutant
lineages are sampled are presumably very rare (e.g.,
Hartl et al. 1994). Yet despite their rarity, these sites
have a large impact on maximum-likelihood estimationof g. Because g enters the likelihood function in thefactor e 2gx, changing g has a much larger impact on thelikelihood of sites with many mutant nucleotides thanthose with few. In other words, a single site with a highfrequency of mutant nucleotides is very strong evidence
for positive selection or low jgj, whereas a site with a low
frequency of mutant nucleotides is not very strong
evidence for the opposite. Thus even a few sites atPolymorphism Frequency Spectrum

which multiple mutant lineages are sampled can cause
large inaccuracies in the inferred g when they areincorrectly assumed by the original PRF method to be
the result of a single high-frequency mutant lineage.
These inaccuracies in g then force corresponding
inaccuracies in the inferred us.
Previous simulation studies have not observed these
problems with the PRF method (Bustamante et al.
2001). However, such simulations themselves implicitly
assumed infinite sites, and hence they cannot be used totest this aspect of the PRF method (to be fair, thesesimulations were done with typical eukaryotic populations in mind, where these problems are much lesssevere). In this article, by contrast, we have simulated afinite number of sites that evolve according to the
Wright–Fisher model.
The role of linkage: While our focus has been on theinfinite-sites approximation, both the original PRF andall of our methods make another key assumption: thateach site is independent of all the others (i.e., freerecombination). This is essentially an assumption of
linkage equilibrium between all polymorphic sites. Thispotentially crucial assumption is likely to be violated inmany real populations, particularly when the method isapplied to estimate selective pressure on short stretchesof DNA (e.g., a single gene) that are linked over longtimescales. This is problematic for both the original PRF
and our methods, but it may be particularly relevant in
populations in which finite-sites issues are important,because recombination is presumably less common inbacteria than in eukaryotes and because higher us
implies a higher density of polymorphic sites, and hencetighter linkage. Ideally, we would address this assumption using theory that allows us to infer the strength of
selection acting on a large number of sites with anarbitrary degree of linkage. However, this is a formidable challenge, and no such theory yet exists.
Despite this, models that assume free recombination
are still useful for several reasons. First, when selection
pressures are weak, sites segregate over long timescales,
so recombination may be frequent enough even among
intragenic sites (or in bacterial genomes) that segregating sites are unlinked over these timescales. Since both
PRF and diffusion-type methods are primarily usefulat estimating selection pressures when they are of order
1/N, the selection pressures our methods can resolve
get smaller as us increases. Thus while populations withlarge us may tend to have lower recombination rates
among polymorphic sites, the polymorphic sites we aresensitive to segregate over longer timescales, and recombination has longer to act. Whether this means thatlinkage can be neglected in any particular population isan empirical question, but fortunately it is usuallystraightforward to estimate linkage disequilibria directlyfrom the data. Thus in any particular case we canestimate whether or not the assumption of free recombination is valid for the purposes of the original PRF

2187

or our modified methods, before applying either
approach. When there is no linkage disequilibrium inthe data, we can use our methods with confidence.
When on the other hand it appears from the data that
free recombination is not a reasonable assumption, amodel that assumes no linkage may still be useful as anull model and a limiting case, and it can be comparedwith more complicated possibilities.
Finally, it is important to note that linkage between
segregating sites will not bias estimates of the selection
pressure g, provided g is equal across sites (Akashi and
Schaeffer 1997; Bustamante et al. 2001). Rather, theprimary effect of linkage is to increase the variance insuch estimates above the predictions of our method(Bustamante et al. 2001). Thus while the predictions ofour models (or the original PRF) will tend to have smaller
confidence intervals than they should, the estimatedvalues of g and us should not be systematically biased. It is
important to note, however, that when linkage is important we must be cautious about testing hypotheses
using these methods; underestimates of the confidence
intervals could, for example, lead us to erroneously reject a neutral null. Further, when selection pressures varysignificantly across sites, the mean parameter estimatesthemselves could be biased; this more complex situation
remains an important topic for future exploration.
Comparison between the PRF and diffusion methods: Our per-site version of the PRF model relaxes some
aspects of the infinite-sites approximation, but it still assumes the mutant lineages do not interact. The diffusion
approach removes the infinite-sites approximation altogether. Thus, the diffusion approach contains none of
the biases associated with infinite-sites approximation
associated with the traditional PRF and, to a lesser
degree, the per-site PRF.
The diffusion method is also easily extendable tomore complex evolutionary situations. For example, we
can explore different selective costs for different nucleotides, or more than one preferred nucleotide. Thesepossibilities lead to obvious modifications of the diffusion equations and their solutions and hence to the
maximum-likelihood estimation. Mutational biases are
more difficult to explore, as no simple steady-state solution to the three-dimensional diffusion equation exists
for biased mutation rates, but computational results canstill be obtained. It is also straightforward to investigatebalancing selection or the effects of dominance: theselead to well-understood modifications to the diffusion
equations and their steady-state solutions (Ewens 2004).
In the PRF framework, by contrast, such generalizations
are much more complex. In particular, balancing selection is impossible to analyze within the PRF framework,
because it leads to mutant lineages reaching stableintermediate frequencies in the population. As a result,the generation of new mutations is not balanced by theextinction or fixation of older ones, and hence nosteady-state distribution of lineage frequencies exists.

2188

M. M. Desai and J. B. Plotkin

The methods we have developed also allow us to relax
the assumption of constant g across sites and instead
infer aspects of the distribution of selection pressures(Piganeau and Eyre-Walker 2003, Nielsen et al. 2005,
and Boyko et al. 2008 have recently analyzed similar
extensions of the original PRF method). It is not yet
clear how much data is required to provide adequate
power for inferring this distribution to a given resolution. However, we can hope to gain a great deal of
insight with only a few additional parameters—say, thenumber of sites that are neutral, lethal, and negatively
selected and the weighted average selection pressure on
the latter class. As we have shown, we can estimate thenumber of lethal sites with no reduction in power
relative to the original PRF method, so this proposal
would involve only one additional parameter. Additional classes of sites would involve more parameters.
Eyre-Walker and Keightley (2007) have recently
stressed that the distribution of fitness effects ofdeleterious mutations is likely to be complex andmultimodal. Using our approach, we can choose how
to focus the power in the data to investigate the aspectsof this complex distribution that are most interesting ina given situation. The appropriate choice of resolution
will depend on the context and quality of the data.
However, the diffusion approach is not without drawbacks. The one-dimensional version does not make useof all of the data available in the observed polymorphism spectrum and is likely to be of little use in practicebecause we typically do not know whether or not
the ancestral nucleotide was preferred. The threedimensional version uses all the data in the folded case,but cannot make use of outgroup data (the unfoldedcase). In practice, when we do not have an outgroup, thethree-dimensional diffusion method provides the mostaccurate estimates of selection pressures across the fullrange of parameters. When an outgroup is available, thethree-dimensional diffusion approach is still sometimes
more accurate than the PRF, despite wasting information on the ancestral state, but only in parameter
regimes typical of bacteria or viruses. In other cases
the per-site PRF method is best.
The diffusion method also cannot naturally handle
positive selection. The steady-state evolutionary dynamics at a site are always dominated by the preferred
nucleotide (or nucleotides), with negative selection
acting against polymorphisms for the disfavored nucleotides. At the level of an individual site, positive selection
is a process that is intrinsically out of steady state: thespread of a favorable nucleotide before it becomes
fixed. The PRF method handles positive selection by
implicitly positing rather strange dynamics at individualsites. All mutant lineages are assumed to be positively
selected—so if a mutant nucleotide fixes at a given site,
mutations back to the ancestral nucleotide are again
assumed to be positively selected (the analogous strange
dynamics also apply to negative selection). While thisassumption makes little sense at a per-site level, it allows
the PRF model to obtain a steady state across sites,
provided positive selection is ongoing and not saturated. We could modify our diffusion methods to mimic
the PRF treatment of positive selection by changing theboundary conditions in our diffusion equations. Specifically, we would assume that probability flowing into x ¼1 (i.e., fixation of a mutant nucleotide) is absorbed andmoved to x ¼ 0 (i.e., ‘‘reset’’ so that new mutations willagain be favored). This diffusion equation can be solved
exactly, and the solution used as a basis for inferring
positive selection using the the per-site diffusion methods we have developed.
Ideally, however, we want to infer positive selection in
the context of a realistic and well-defined model of the
dynamics at individual sites. Such an approach would
necessarily involve full time-dependent solutions to the
diffusion equations; Mustonen and Lassig (2007)
suggest a method along these lines. We do not pursue
this approach here, but see Chen et al. (2007) for a step inthis direction. Regardless of the methodology, it will
always be difficult to discriminate positive selection fromnegative selection on the basis of the polymorphism
frequency spectrum alone, particularly when only folded
data are available. Whether selection is positive or
negative, mutant lineages drift nearly neutrally when
their frequency is between 0 and 1/jgj. Positively selected
lineages then fix relatively quickly once their frequency
becomes substantially larger than 1/jgj, while negatively
selected lineages rarely ever reach frequencies larger
than 1/jgj. Thus from the point of view of the polymorphism frequency spectrum, positive selection issimilar to random drift on [0, 1/jgj], with the upperbound a roughly absorbing boundary condition. Negative selection, on the other hand, is also similar to random
drift on [0, 1/jgj], but with the upper bound a roughlyreflecting boundary condition. Although this is relatively
crude—selection does in fact have some impact on lowfrequency lineages, and the boundary conditions are notexactly absorbing or reflecting—it indicates that thepolymorphism frequency spectrum is roughly similar
for negative and positive selection at the same jgj. Thus
power to distinguish positive from negative selection
based on the polymorphism frequency spectrum, especially with folded data, will always be relatively limited,regardless of the method used.
Different steady states: The PRF framework and the
diffusion methods are based on fundamentally different
steady states. The PRF assumes an ancestral state; thesteady-state frequency distribution is that of derived
nucleotides relative to this ancestral state. By contrast,
the diffusion methods assume mutation back and forth
between all possible nucleotides at each site. The steadystate frequency distribution is a full mutation-selectiondrift balance of all four nucleotides at the site. It has‘‘forgotten’’ the ancestral state, and it is obtained over alonger timescale than the PRF steady state.

Polymorphism Frequency Spectrum

The differences in the nature of the two steady states
are most readily apparent when considering their
behavior for small us. This is easiest to see by comparing
the PRF steady state to the one-dimensional diffusion
result (assuming as usual that one nucleotide is preferred and the other three are unpreferred). For small
us, the population will typically be fixed for one of thefour nucleotides. Occasionally mutations will occur,
making the site temporarily polymorphic within the
population. Usually these mutant lineages will drift at
low frequency for a short while before going extinct, but
occasionally they will fix and the population will be fixed
for this new nucleotide. The PRF steady state describes
the (transient) frequency spectrum of these occasional
mutant lineages before they fix or go extinct. This
spectrum is always measured relative to the ancestral
nucleotide, which is the last mutation that fixed. On
the other hand, the diffusion method describes thelonger-term steady-state frequency distribution of thenucleotides. This includes some weight near x ¼ 0(the preferred nucleotide is fixed, or nearly so), as well
as some weight near x ¼ 1 (one of the unpreferred
nucleotides is fixed, or nearly so). Because of theselective bias, the population will be fixed or nearly
fixed for the preferred nucleotide a fraction 1=ð1 1 3e 2g Þof the time (Ewens 2004), so 1=ð1 1 3e 2g Þ of the total
weight is in the part of the distribution near x ¼ 0. The
remainder of the time, the population will be fixed ornearly fixed for one of the unpreferred nucleotides; thisis the weight in the distribution near x ¼ 1.
The discussion above helps to elucidate the relationship between the steady states in the PRF model vs. thosein the diffusion model. A fraction 1=ð1 1 3e 2g Þ of thetime the most recent fixation will have been for the
preferred state. This is the ancestral state for the PRF
method, and the frequency spectrum of the occasional
unpreferred negatively selected mutants is given by the
PRF steady state. The remainder of the time the mostrecent fixation will have been for one of the unpreferred states. Now this is the ancestral state for the PRF
method, and the frequency spectrum of occasional
mutants for the preferred nucleotide, which are underpositive selection, is again given by the PRF steady state,
but with the opposite sign of g. Thus we expect that the
diffusion method steady state is just the sum of these two
PRF steady states, weighted by the frequency with which
each type of nucleotide is the ancestral state, namely
f diff ðg; xÞ  pðgÞf PRF ðg; xÞ 1 ½1  pðgÞf PRF ðg; 1  xÞ;

where p(g) denotes the chance that the most recent
fixation was for the preferred state. This relationshipbetween the two steady states is illustrated in Figure 3. It
holds in the limit of small us. For larger us, thepopulation is not always nearly fixed for one state.
Mutant lineages occur more frequently and can segregate simultaneously (i.e., finite-sites effects). Since thediffusion approach accounts for this while the PRF does

2189

Figure 3.—Relationship between the PRF and diffusion
frequency spectra for small us. Shown at the top is the one-dimensional diffusion frequency spectrum for us ¼ 103, g ¼
0.5. Most of the time the population is fixed or nearly fixed
for the preferred nucleotide (this is the weight near x ¼ 0),
while sometimes the population is fixed or nearly fixed for
an unpreferred nucleotide (the weight near x ¼ 1). Whenthe population is nearly fixed for the preferred (unpreferred)
nucleotide, mutations produce a PRF site frequency spectrum(SFS) as shown on the left (right) of the bottom graph. For
this small value of us, the diffusion frequency spectrum isroughly equal to the sum of these PRF frequency spectra.

not, the relationship between the steady states breaks
down.
The difference between the diffusion and PRF steadystates points to one additional bias in the PRF methodthat we have not yet discussed. Imagine a stretch of
sequence that has been under constant purifying
selection for an extremely long time. As we have justseen, we expect that at a fraction 3e 2g =ð1 1 3e 2g Þ of thesites in this sequence, the most recent ancestor will be
an unpreferred state. The PRF frequency spectrum atthese sites will be characteristic of positive selection(mutations to the preferred nucleotides are beneficial).
Although this will typically represent a minority of thesites, the PRF method weights sites under positive
selection more strongly than those under negative
selection when inferring a weighted average g. Thus,
the PRF method will be biased toward inferring positive
selection—despite the fact the sequence has experienced negative selection for an arbitrarily long time.
This effect is present even for small us, and it is most
severe when selection is weak but nonzero.
Whether this effect represents a bias in the PRF
method is a matter of taste. In some sense, the PRF is2190

M. M. Desai and J. B. Plotkincorrectly identifying sites under positive selection: after
all, the mutations to the preferred nucleotides are infact beneficial. In another sense, however, the sequence
has experienced constant purifying selection throughout its past. The fact that occasionally an unpreferred
nucleotide fixes before reverting to the preferred state
does not reflect any change in selection pressures. Thus,
if we define positive selection as a pressure to innovate,
or a shift in the selective landscape, then such areversion does not reflect positive selection. We have
explored this effect through simulation and found thatin microbial parameter regimes it is a detectable but
minor effect and is much less important in practice thanthe other biases we have studied. However, it would bestraightforward to modify either the original or the persite PRF method to correct for this (the diffusion
methods already account for this effect). We simply
change the expected spectrum to reflect the fact that afraction 1=ð1 1 3e 2g Þ of sites will exhibit the spectrum
characteristic of negative selection while a fraction
3e 2g =ð1 1 3e 2g Þ will exhibit that characteristic of positive
selection (with equal but opposite g). This leads to asimple modification of the PRF spectrum, with noadded parameters, which we can use for ML estimation
as before.
Different timescales: Because the PRF frequency
spectrum is measured relative to the last fixed state, itreaches steady state in at most of order N generations,the coalescent timescale. The diffusion method takes
longer to reach steady state—in the neutral case,mutations fix every 1/m generations, so it takes of order
1=m generations to reach steady state. When selection is
operating, the diffusion steady state can arise either
more quickly or more slowly, depending on the initialstate of the population. Note that the ratio of the PRF to
the diffusion timescale is of order us. Thus for the small
us relevant for most eukaryotic populations, there is a
big gap between the two timescales, while for the larger
us relevant to bacteria and viruses the two timescales will
be more similar.
These two timescales determine whether the presence of information from an outgroup is likely to
improve our inferential power. Imagine we sample
polymorphism from some population A and one sequence from an outgroup species. If the outgroupdiverged from the population A less than N generations
ago, population A has not reached the PRF steady statein the time since the divergence. Thus the frequency
spectrum of population A is not independent from thefrequency spectrum in the outgroup. The outgroupdoes not provide much more information than simply
sampling another individual from population A(though it does provide some additional information).
This situation has recently been analyzed in the neutralcase by Chen et al. (2007) and applied to humanpolymorphism data using a Neanderthal outgroup.

If instead the outgroup diverged from population
A more than N but less than 1/m generations ago, it canbe used to infer the most recent ancestral state. Since
the time of this ancestor, the PRF frequency spectrum
will have reached equilibrium, and it can be used toestimate parameters accounting for this knowledge(though finite-sites issues can still matter when us islarge, and diffusion methods may be more accurate
despite the information from the outgroup). But thecloser the time since the outgroup diverged is to 1/m
generations, the more unreliable the inference of theancestral state becomes, and ancestral state misidentification becomes more of a problem. This problem andpossible ways to account for it have recently been
examined by Hernandez et al. (2007).
If the outgroup diverged more than 1/m generations
ago, the population will have forgotten the ancestral
state since divergence, and the outgroup will not be
useful. Of course, in this case we would probably noteven be able to align the outgroup sequence to thesequences from population A. In this case, only thefolded frequency data are available, and the diffusion
methods are most appropriate.
We thank Daniel Fisher, Marc Feldman, and Warren Ewens forhelpful comments on the manuscript. M.M.D. acknowledges support
from Center grant P50GM071508 from the National Institute ofGeneral Medical Science to the Lewis–Sigler Institute. J.B.P. acknowledges support from the Burroughs Wellcome Fund, the James S.
McDonnell Foundation, and the Defense Advanced Research Projects
Agency Fun-Bio Program (HR0011-05-1-0057).

LITERATURE CITED
Akashi, H., 1995 Inferring weak selection from patterns of polymorphism and divergence at ‘‘silent’’ sites in Drosophila DNA. Genetics 139: 1067–1076.
Akashi, H., 1999 Inferring the fitness effects of DNA mutations
from polymorphism and divergence data: statistical power to detect directional selection under stationarity and free recombination. Genetics 151: 221–238.
Akashi, H., and S. W. Schaeffer, 1997 Natural selection and thefrequency distributions of ‘‘silent’’ DNA polymorphism in Drosophila. Genetics 146: 295–307.
Bartolome, C., X. Maside, S. Yi, A. L. Grant and B. Charlesworth,2005 Patterns of selection on synonymous and nonsynonymous
variants in Drosophila miranda. Genetics 169: 1495–1507.
Boyko, A. R., S. H. Williamson, A. R. Indap, J. D. Degenhardt, R. D.
Hernandez et al., 2008 Assessing the evolutionary impact ofamino acid mutations in the human genome. PLoS Genet. 4:e1000083.
Bustamante, C. D., J. Wakeley, S. Sawyer and D. L. Hartl,
2001 Directional selection and the site-frequency spectrum. Genetics 159: 1779–1788.
Bustamante, C. D., R. Nielsen, S. Sawyer, K. M. Olsen, M. D.
Purugganan et al., 2002 The cost of inbreeding in Arabidopsis.
Nature 416: 531–534.
Bustamante, C. D., R. Nielsen and D. L. Hartl, 2003 Maximum
likelihood and Bayesian methods for estimating the distributionof selective effects among classes of mutations using DNA polymorphism data. Theor. Popul. Biol. 63: 91–103.
Bustamante, C. D., A. Fledel-Alon, S. Williamson, R. Nielsen,
M. T. Hubisz et al., 2005 Natural selection on protein-coding
genes in the human genome. Nature 437: 1153–1157.
Chen, H., R. E. Green, S. Paabo and M. Slatkin, 2007 The joint allelefrequency spectrum in closely related species. Genetics 177: 387–398.

Polymorphism Frequency Spectrum
Desai, M. M., and D. S. Fisher, 2007 Beneficial mutation-selection
balance and the effect of linkage on positive selection. Genetics
176: 1759–1798.
Drake, J. W., B. Charlesworth, W. Charlesworth and J. F. Crow,
1998 Rates of spontaneous mutation. Genetics 148: 1667–
1686.
Ewens, W. J., 2004 Mathematical Population Genetics: I. Theoretical Introduction. Springer, New York.
Eyre-Walker, A., and P. D. Keightley, 2007 The distribution of fitness effects of new mutations. Nat. Rev. Genet. 8: 610–618.
Galtier, N., E. Bazin and N. Bierne, 2006 Gc-biased segregation of
noncoding polymorphisms in Drosophila. Genetics 172: 221–
228.
Hartl, D. L., E. N. Moriyama and S. A. Sawyer, 1994 Selection intensity for codon bias. Genetics 138: 227–234.
Hernandez, R. D., S. H. Williamson and C. D. Bustamante,2007 Context dependence, ancestral misidentification, andspurious signatures of natural selection. Mol. Biol. Evol. 24:
1792–1800.
Lercher, M. J., N. G. C. Smith, A. Eyre-Walker and L. Hurst,
2002 The evolution of isochores: evidence from SNP frequency
distributions. Genetics 162: 1805–1810.
Loewe, L., B. Charlesworth, C. Bartolome and V. Noel,2006 Estimating selection on nonsynonymous mutations. Genetics 172: 1079–1092.
Lynch, M., and J. S. Conery, 2003 The origins of genome complexity. Science 21: 1401–1404.
McDonald, J., and M. Kreitman, 1991 Adaptive protein evolution
at the adh locus in Drosophila. Nature 351: 652–654.
Moran, P. A. P., 1959 The survival of a mutant gene under selection.
ii. J. Aust. Math. Soc. 1: 485–491.
Mustonen, V., and M. Lassig, 2007 Adaptations to fluctuating selection in Drosophila. Proc. Natl. Acad. Sci. USA 104: 2277–2282.
Nachman, M. W., 1998 Deleterious mutations in animal mitochondrial DNA. Genetica 102/103: 61–69.
Nielsen, R., C. D. Bustamante, A. G. Clark, S. Glanowski, T. B.
Sackton et al., 2005 A scan for positively selected genes in
the genomes of humans and chimpanzees. PLoS Biol. 3: 976–
985.
Piganeau, G., and A. Eyre-Walker, 2003 Estimating the distribution of fitness effects from DNA sequence data: implications
for the molecular clock. Proc. Natl. Acad. Sci. USA 100:
10335–10340.
Rand, D. M., and L. M. Kann, 1998 Mutation and selection at silent
and replacement sites in the evolution of animal mitochondrialDNA. Genetica 102/103: 393–407.
Sawyer, S., R. J. Kulathinal, C. D. Bustamante and D. L. Hartl,
2003 Bayesian analysis suggests that most amino acid replacements in Drosophila are driven by positive selection. J. Mol. Evol.
57: S154–S164.
Sawyer, S. A., and D. L. Hartl, 1992 Population-genetics of polymorphism and divergence. Genetics 132: 1161–1176.
Wakeley, J., 2003 Polymorphism and divergence for island-model
species. Genetics 163: 411–420.
Watterson, G. A., 1977 Heterosis or neutrality? Genetics 85: 789–
814.
Weinreich, D. M., and D. M. Rand, 2000 Contrasting patterns ofnonneutral evolution in proteins encoded in nuclear and mitochondrial genomes. Genetics 156: 385–399.
Williamson, S., A. Fledel-Alon and C. D. Bustamante,
2004 Population genetics of polymorphism and divergencefor diploid selection models with arbitrary dominance. Genetics
168: 468–475.
Williamson, S. H., R. Hernandez, A. Fledel-Alon, L. Zhu, R. Nielsen
et al., 2005 Simultaneous inference of selection and population
growth from patterns of variation in the human genome. Proc.
Natl. Acad. Sci. USA 102: 7882–7887.
Wright, S., 1938 The distribution of gene frequencies underirreversible mutation. Proc. Natl. Acad. Sci. USA 24: 253–
259.

2191

Wright, S., 1949 Adaptation and selection, pp. 365–389 in Genetics,
Paleontology, and Evolution, edited by G. L. Jepson, G. G. Simpson
and E. Mayr. Princeton University Press, Princeton, NJ.
Zhu, L., and C. D. Bustamante, 2005 A composite-likelihood approach for detecting directional selection from DNA sequence
data. Genetics 170: 1411–1421.
Communicating editor: M. K. Uyenoyama

APPENDIX: ALGORITHM TO COMPUTE THE
LIKELIHOOD EXPRESSIONS UNDER
THE PER-SITE PRF

To implement an inference procedure based on our
per-site version of the PRF model we need to compute
likelihood expressions for any possible unfolded configuration of nucleotides, an example of which is shownin Equation 10.
Our algorithm starts by initializing the expressions
associated with each configuration to zero. We loop overall partitions of the sample size n written as the sum ofnonzero integers, n ¼ n1 1 n2 1 . . . 1 nz 1 nw. In thispartition, numbers n1 through nz represent the numberof samples of distinct mutant lineages, and nw represents the number of samples with the ancestral state.
The sum of all such sample numbers equals the total
sample size, n. Such a sample will correspond to aparticular folded configuration (a, b, c, d) that depends
upon the identity of the mutant nucleotides associatedwith the z mutant lineages. There are 3z distinct possible
assignments of mutant nucleotide types to the mutant
lineages. We loop over each of these possible assignments, and for each we add a term of the form
ð1=3ÞzF ð1Þn1 F ð2Þn2
F ðnz Þnz
3

n1 !
n2 !
nz !

to the likelihood expression for the folded configuration (a, b, c, d) that arises from the assignment of
nucleotide types. After looping over all possible numbers of mutant lineages and over all possible assignments of mutant nucleotide types to such lineages, we
obtain complete expressions for the likelihoods of all
unfolded configurations.
To compute likelihood expressions for folded configurations, we perform the exact same algorithm as
above with a small modification. In this case, once wehave assigned nucleotide types to each of the z mutant
lineages, we determine the associated folded configuration by simply sorting the number of lineages (including the wild type) of each mutant nucleotide type indescending order.
A computer implementation of this algorithm, foreither folded or unfolded likelihood expressions, isfreely available upon request.