New Explicit Expressions for Relative Frequencies of Single-Nucleotide Polymorphisms With Application to Statistical Inference on Population Growth
A. Polanski, M. Kimmel

Abstract

We present new methodology for calculating sampling distributions of single-nucleotide polymorphism (SNP) frequencies in populations with time-varying size. Our approach is based on deriving analytical expressions for frequencies of SNPs. Analytical expressions allow for computations that are faster and more accurate than Monte Carlo simulations. In contrast to other articles showing analytical formulas for frequencies of SNPs, we derive expressions that contain coefficients that do not explode when the genealogy size increases. We also provide analytical formulas to describe the way in which the ascertainment procedure modifies SNP distributions. Using our methods, we study the power to test the hypothesis of exponential population expansion vs. the hypothesis of evolution with constant population size. We also analyze some of the available SNP data and we compare our results of demographic parameters estimation to those obtained in previous studies in population genetics. The analyzed data seem consistent with the hypothesis of past population growth of modern humans. The analysis of the data also shows a very strong sensitivity of estimated demographic parameters to changes of the model of the ascertainment procedure.

A lot of research has been done to develop methods for discovery of single-nucleotide polymorphisms (SNP) and to characterize distributions of SNPs across the genome (Collinset al. 1997; Wanget al. 1998; Cargillet al. 1999; Marthet al. 1999; Picoult-Newberget al. 1999; Altshuleret al. 2000). SNP data have already been used in association studies of complex diseases (Boerwinkleet al. 1996; Halushkaet al. 1999; Bonnenet al. 2000; Trikkaet al. 2002), and it is believed that eventually they will enable creation of fine genetic maps for complex traits analysis (Kruglyak 1999; Rish 2000). Databases, like that of the SNP Consortium, at http://snp.cshl.org, contain massive amounts of data on positions of SNPs in the human genome, but it is likely that most of these SNPs are very rare and therefore of limited value in gene association studies. Estimates of distributions of expected relative frequencies of SNPs result from studies that use population genetics models, e.g., Durrett and Limic (2001) and Wang et al. (1998), and the predicted excess of rare alleles is explained as resulting from expansion of human populations.

Using population genetics methods to model and analyze SNPs opens an area for investigating problems like predicting frequencies of SNPs under various demographic scenarios, inferring demographic parameters and history from sampling frequencies of SNPs, comparing estimates obtained on the basis of SNP data to those obtained with other methods, and evaluating efficiency of using SNP data for estimation of population parameters. Several interesting studies were carried out in this area. Studies by Durrett and Limic (2001) and Wang et al. (1998) estimated frequencies of SNPs under the hypothesis of population growth. A problem of how sampling frequencies of SNPs are influenced by ascertainment procedures was investigated by Eberle and Kruglyak (2000), Yang et al. (2000), and Renwick et al. (2002). Using SNPs for estimation of the scaled product parameter θ= 4Neμ of the effective population size Ne and mutation rate μ, under assumption of constant population size, was studied by Kuhner et al. (2000). They took into account various hypotheses of spatial (chromosomal) distributions of SNPs such as complete or partial linkage or occurrence of linked segments of nonrecombining SNPs and, on the basis of extensive simulations, evaluated accuracy of estimates and possible sources of bias. Studies by Nielsen (2000) and Wakeley et al. (2001) were devoted to detection of signatures of human population growth in SNP data. Nielsen (2000) fitted the scenario of exponential expansion to SNP data of Picoult-Newberg et al. (1999). Wakeley et al. (2001) used the model of stepwise change of the population size with population subdivision (Wakeley 2001). They fitted their model to SNP data from Wang et al. (1998), Cargill et al. (1999) and Altshuler et al. (2000). Parameter-space regions corresponding to the highest likelihoods were not inconsistent with the hypothesis of population growth. Moreover, if the ascertainment bias was not considered, less realistic shapes of parameter regions were obtained. Comparison of cases in which population substructure was not considered with those in which it was considered seems to support the latter scenario. To evaluate SNP frequencies, these studies used the standard coalescence approach and Monte Carlo simulations.

Sampling distributions of SNP frequencies in populations with time-varying size can be calculated with the use of analytical expressions for the expected lengths of branches in the coalescence tree derived in the articles by Griffiths and Tavare (1998), Wooding and Rogers (2002), and Polanski et al. (2003). Analytical expressions allow computations, which are faster and more accurate than Monte Carlo simulations. However, the approaches shown in the articles by Griffiths and Tavare (1998), Wooding and Rogers (2002), and Polanski et al. (2003) suffer from one common difficulty, numerical instability for larger genealogies. When the analyzed genealogy size is >50, these analytical methods are either inapplicable or difficult to apply, due to the explosion of coefficients in equations. Wooding and Rogers (2002) give a method to stabilize numerical computations, which is valid for the case where effective population size changes in a stepwise manner. Here we show another approach, which is more general in the sense that it does not require assumption of piecewise constant history of effective population size. We transform equations for the relative frequencies of SNP to the form with nondiverging coefficients and we provide expressions, obtained with the use of methods of hypergeometric summation, to compute these coefficients. We also provide analytical expressions to describe the influence of the discovery procedure (ascertainment) on SNP frequencies. Our methods allow us to perform tasks that otherwise are prohibitive or cumbersome, like analyzing large genealogies, estimating confidence limits for parameters by resampling studies, and studying sensitivity of models to parameter changes. Using our methods we study our power to test the null hypothesis of evolution with constant population size vs. the alternative hypothesis of population expansion, for SNP data, under the exponential model of population size change. We also analyze some of the available SNP data (Picoult-Newberget al. 1999; Trikkaet al. 2002) and we compare our results to those obtained in previous studies concerning estimation of populations size changes (Slatkin and Hudson 1991; Rogers and Harpending 1992; Polanskiet al. 1998; Weiss and Haeseler 1998).

METHODS

We consider the process of coalescence with time-changing effective population size. Notation for the coalescence tree, for the sample of size n = 5 DNA sequences, is shown in Figure 1. Time t is measured, in number of generations, from the present to the past. Random times between coalescence events are denoted by Sn, Sn–1,..., S2 and sn, sn–1,..., s2. Cumulative times to coalescence, from sample of size n to sample of size k – 1, are denoted by Tk, k = 2, 3,..., n, and their realizations by corresponding lowercase letters tn, tn–1,..., t2, 0 < tn < tn–1... < t2.

Figure 1.

—Notation for ancestral history of a sample of DNA sequences. Coalescence times for the sample of size n = 5 are denoted by T5, T4,..., T2 and their realizations by corresponding lower case letters t5, t4,..., t2. Times between coalescence events are denoted by S5, S4,..., S2 and s5, s4,..., s2. A mutation event is marked by an open circle. Sequences 4 and 5 have mutant alleles (bases), while sequences 1–3 have ancestral ones.

We assume that an observed SNP was produced by a single, neutral mutation, like the one denoted in Figure 1 by an open circle. In Figure 1 sequences 4 and 5 have mutant alleles (bases), while sequences 1, 2, and 3 have ancestral ones. In the situation where it is not known which allele is mutant and which is ancestral, we use the terms rare and frequent allele. In other words, the SNP in Figure 1 has configuration b = 2 mutant vs. nb = 3 ancestral, or b = 2 rare vs. nb = 3 frequent alleles. We assume that mutation intensity for SNPs is very low; i.e., they follow the infinite-sites mutation model.

Probability that a SNP has b mutant bases: Probability qnb that a SNP site in a sample of n chromosomes has b mutant bases, under the infinite-sites mutation model, is given by Griffiths and Tavare (1998, Equation 1.3) in terms of expectations of times in the coalescence tree (see also articles by Fu 1995; Sherryet al. 1997; Nielsen 2000; Wooding and Rogers 2002). In our notation, this expression has the form qnb=((nb1)!(b1)!(n1)!)k=2nk(k1)(nkb1)E(Sk)k=2nkE(Sk), (1) where 0 < b < n, Sk = TkTk+1, and Tn+1 = 0.

The above expression can be written as qnb=((nb1)!(b1)!(n1)!)k=2nj=knj(j1)(nkb1)Akjnejk=2nj=knj(j1)(k1)Akjnej (2) (Polanskiet al. 2003), where ej=0tqj(t)dt (3) are expectations of times distributed as qj(t)=(j2)Ne(t)exp(0t(j2)dσNe(σ)), (4) with the effective population size history described by a function of reverse time, Ne(t),t[0,). (5) Coefficients Akjn are defined by the expression Akjn=l=k,ljn(l2)l=k,ljn[(l2)(j2)],kjn, (6) Annn=1 .

Equation 2 is an analytic expression for probabilities qnb. Wooding and Rogers (2002) derive equations with the structure analogous to Equations 2, 3, 4, 5, which also contain expectations defined in (3). In contrast to (2), they do not provide explicit expressions for coefficients in the equations; instead they use linear algebra operations (matrix diagonalization) to compute them. Both articles (Wooding and Rogers 2002; Polanskiet al. 2003) report that it is rather difficult to efficiently apply analytical formulas for genealogies of size n > 50 because of the occurrence of diverging terms with alternating signs.

Methods for computation of qnb for large genealogies: To avoid large numerical errors in summations in (2) for genealogies n > 50, one needs to apply computations with precision of hundreds, or even thousands, of decimal digits (Wooding and Rogers 2002), which significantly slows down computational process and requires appropriate software. Such computations must be also carefully executed. It is necessary to repeat computations several times, with an increasing accuracy, and to examine the convergence of the returned values.

Wooding and Rogers (2002) developed a way to avoid the necessity of extending precision of the arithmetics, based on a uniformization technique of computing matrix exponents. It is applicable for the case when the population size changes in a stepwise (piecewise constant) manner, with a finite number of steps, and it allows evaluating the expressions in a standard double precision arithmetic. However, when the number of steps in the population size history becomes large, e.g., if one attempts to approximate a given Ne(t) by a piecewise constant function, this approach may be difficult to apply.

Below, we present a method for computing qnb for large genealogies, which is more general than the one developed by Wooding and Rogers (2002), in the sense that it does not require assumption of stepwise change of Ne(t). The idea is to reverse the order of summation in both denominator and numerator in Equation 2. We observe that the resulting expressions contain coefficients that do not explode when n increases. Proceeding in this way we obtain qnb=j=2nejk=2jj(j1)(nkb1)((nb1)!(b1)!(n1)!)Akjnj=2nejk=2jj(j1)(Akjn(k1)) (7) =j=2nejWbjnj=2nejVjn. (8) In the above, we introduced coefficients Vjn=k=2jj(j1)Akjnk1 (9) and Wbjnk=2jj(j1)(nkb1)(nb1)!(b1)!(n1)!Akjn. (10) For k > nb + 1, the elements in the sum (10) become zero, so the upper limit j can be replaced by min(j, nb + 1). Coefficients Vjn and Wbjn remain the same for all histories of effective population size Ne(t). Once calculated, they can be stored in computer memory or tabularized and reused when we wish to analyze different histories Ne(t), e.g., when maximizing likelihood function with respect to population growth parameters. Their important property is that their growth, when genealogy size n increases, is very moderate; e.g., for n = 100, maxj(Vj100)=17.13 , maxb,jWbj100=8.24 ; for n = 500, maxj(Vj500)=38.36 , maxb,jWbj500=18.36 ; and for n = 1000, maxj(Vj1000)=54.18 , maxb,jWbj1000=25.94 . In Figure 2 growth plots of maxb,jWbjn and maxj(Vjn) vs. n are shown. One can see that both plots are, asymptotically, of the power type with the exponent less than one.

Figure 2.

—Growth plots of Formula and Formula (○) vs. n.

Expressions in (10) and (9) are sums of hypergeometric series, which can be seen by factoring the denominators in (6), (2l)(2j)=12(lj)(l+j1) , and then expressing coefficients Akjn in (6) as follows: Akjn=n!(n1)!(n+j1)!(nj)!(2j1)j(j1)(j+k2)!(k1)!(k2)!(k2)!(jk)!(1)jk. (11) Substituting (11) in (9) and using Chu-Vandermode identity (Grahamet al. 1998, p. 212, Equation 5.93) we obtain Vjn=(2j1)n!(n1)!(n+j1)!(nj)![1+(1)j]. (12) Coefficients Wbjn in expression (10) can be efficiently computed with the use of recursive procedures (Paule and Schorn 1994; Petkovseket al. 1996). Several recursions for Wbjn are possible, depending on which index one decides to consider as the running one. We used the implementation of Zeilberger's algorithm in Mathematica, developed by P. Paule and M. Schorn, available at http://www.risc.uni-linz.ac.at/research/combinat/risc/software/, to obtain recursions for Wbjn . We found the following recursive scheme, with respect to the index j, very useful: Wb2n=6(n+1), (13) Wb3n=30(n2b)(n+1)(n+2), (14) Wb,j+2n=(1+j)(3+2j)(nj)j(2j1)(n+j+1)Wbjn+(3+2j)(n2b)j(n+j+1)Wb,j+1n. (15) The above recursions are numerically stable and very fast. We used them, implemented in a standard double-precision arithmetic, for genealogies consisting of thousands of DNA sequences (the largest value of n tested was n = 5000). We did not perform precise measurements, but usually, when calculating probabilities qnb, according to (8), computing coefficients Wbjn and Vjn takes only a small fraction of the time, while most of the computing effort is needed to evaluate expectations ej.

Influence of the ascertainment procedure on SNP sampling frequencies: Most of the published data on SNP sampling frequencies are obtained in a two-step process, where the first step involves discovering chromosomal locations of a number of SNPs, and the second one involves DNA sequencing of a sample of n chromosomes restricted to locations discovered in the first step. The first step is called SNP ascertainment and is based on number of chromosomes smaller than n. As described in previous studies, taking into account the ascertainment scheme is a very important aspect of SNP data analysis. Below we derive expressions for modeling the way in which ascertainment modifies SNP sampling frequencies.

We use the following notation introduced by Wakeley et al. (2001): Data set size is n = nD + nO, and ascertainment set size equals to nO + nA, where nA stands for the number of ascertainment-only samples; nO, the number of overlapping samples (both in the ascertainment study and in the later data set); and nD, the number of data-only samples.

To determine how ascertainment modifies probability distribution (22), we merge ascertainment and data sets to obtain the joint set of size nJ = nD + nO + nA. We treat the ascertainment procedure as sampling SNP alleles, without replacement, from the joint set. A SNP is discovered if (a) both alleles are present in the ascertainment sample and (b) none of the alleles in the ascertainment sample has number of copies less than G, where G is a predetermined threshold. Since the joint set contains elements of two types (two alleles), the number of copies of alleles in the ascertainment sample follows a hypergeometric distribution. We analyze two cases: (i) no overlap, which means nO = 0, n = nD, nJ = nD + nA; and (ii) overlap only, which means nA = 0, n = nJ = nD + nO. The case where both overlap and ascertainment-only samples are present is obtained by combining i and ii. We compute frequencies of discovered SNPs in the data set, which follow from conditions a and b above. We analyze first the case i. If a SNP in the joint set has b mutant and nJb ancestral bases, then the probability that a sample of size nA from the joint set has β mutant and nA –β ancestral bases is h(β,nJ,b,nA)=(bβ)(nJbnAβ)(nJnA). (16) For a SNP to be discovered, β must satisfy G ≤β ≤ nAG, with G defined as above. Moreover, the following inequalities must hold: β≤ b, nA –β≤ nJb. Consequently, the probability πnγA that a discovered SNP in the data-only set i has γ= b –β mutant and nD –γ ancestral alleles is πnDγA=β=GnAGqnJγ+βh(β,nJ,γ+β,nA)g=0nDβ=GnAGqnJg+βh(β,nJ,g+β,nA), (17) γ= 0, 1,..., nD. Probabilities qnb are given by (8). The relation γ= b –β follows from the fact that β chromosomes with mutant bases are removed from the joint set. The numerator in (17) is a sum of contributions to πnDγA for possible values of β, while the denominator is a normalizing factor. For case ii assume again that a SNP in the joint set has b mutant and nJb ancestral bases. The probability that a sample of nO has β mutant and nO –β ancestral bases is given by (16) with nA replaced by nO. For this SNP to be discovered β must satisfy G ≤β ≤ nOG. Consequently, the probability πnJbO that a discovered SNP in the joint set ii has b mutant and nJb ancestral alleles is πnJbO=qnJbβ=GnOGh(β,nJ,b,nO)σ GnJGqnJσβ GnOGh(β,nJ,σ,nO), (18) b = G,..., nJG.

If it is not known which of the alleles is mutant and which one is ancestral, we need to symmetrize πnDγA and πnJbO to get probability of data configuration. For case i we have expression P(XR=γ)=pnDγA=πnDγA+πnD,nDγA[1δ(γ,nDγ)],γ=0,1,,[nD2] (19) for the probability that the rare allele has γ copies. For case ii the probability that there are b copies of the rare allele is P(XR=b)=pnJbO=πnJbO+πnJ,nJbO[1δ(b,nJb)],b=G,G+1,,[nJ2]. (20) In the above [n/2] denotes greatest integer ≤n/2.

In the sequel, we refer to the models described above as type i and type ii ascertainment, respectively.

Likelihood function of the sample: Data studied are derived from a number of SNP sites. Let us denote the number of SNP loci by M and random variables defined by diallelic data by [X1,X2,,XM]=[(X1R,X1F),(X2R,X2F),(XmR,XmF),(XMR,XMF)], (21) where XmR is the number of copies of the less frequent (rare) allele and XmF is the number of copies of the more frequent one, in the sample of nm=XmR+XmF . It is possible that XmR=XmF for some indices m, in which case both alleles are equally frequent. We assume that the ancestral state is not known. Then, for an SNP (XmR,XmF) , the probability that we observe configuration bm, nmbm, bm ≤ [nm/2] is P(XmR=bm)=pnmbm=qnmbm+qnm,nmbm[1δ(bm,nmbm)], (22) where δ(·) is the Kronecker delta function and qnb are probabilities defined and evaluated in the previous section.

When SNP sites are located far from one another, random variables {X1, X2,..., XM} in (21) are independent. If the observed numbers of copies of rare alleles are X1R=b1,X2R=b2,,XmR=bm,,XMR=bM , then the log likelihood of the sample (21) is l=m=1Mlog(pnmbm) (23) (Nielsen 2000; Wooding and Rogers 2002). If sample sizes are equal for all SNP loci, n1 = n2 =... = nM = n, the above expression can be written as l=b=1[n2]cblog(pnb), (24) where cb denotes number of SNP loci in the sample, which have configuration of b copies of the rare allele vs. nb copies of the frequent allele. Subsequently, we use expressions (23) and (24) to compute likelihoods of SNP samples with different ascertainment models. To specify the ascertainment model we substitute in (23) or (24), pnb = pnb [expression (22), no ascertainment step], pnb=pnbA [expression (19), ascertainment model type i], or pnb=pnbO [expression (20), ascertainment model type ii].

RESULTS

Exponential history of population size: In our computations we assume an exponential history of effective population size. In previous publications devoted to SNP and demography, Nielsen (2000) assumed an exponential history of Ne(t). However, his analysis is very brief and restricted only to simulations. Others (Wakeleyet al. 2001; Wooding and Rogers 2002) analyzed stepwise histories of effective population size changes.

For an exponential scenario of population growth Ne(t)=Ne0exp(rt), (25) expectations in (3) become ej=ej(Ne0,r)=exp[(j2)(rNe0)1]rEi[(j2)(rNe0)1] (26) (Slatkin and Hudson 1991), where Ei denotes the exponential integral, Ei(μ)=1exp(1μx)xdx,Re(μ)>0 , (Gradshteyn and Ryzhik 1980, Sect. 3.351.5). When the argument (2l)(rNe0)1 in (26) becomes large, computing el(Ne0, r) involves solving product of the type ∞ · 0. For (2l)(rNe0)1>50 , we used expansion, Ei(x)=exp(x)k=1K(1)k(k1)!xk+RK (27) (Gradshteyn and Ryzhik 1980, Sect. 8.215) with RK<K!xK+1, (28) which allowed canceling exp[(2l)(rNe0)1] in (26).

It turns out that sampling frequencies of SNPs depend only on the product parameter κ= rNe0 of initial effective population size and exponential factor.

Distributions of SNP frequencies: Figure 3 provides examples of probabilities of different configurations of SNP sites, for sample size n = 30, and different values of the parameter κ (0, 1, 10), under the assumption that data collection did not include an ascertainment step [expression (22)] or under the ascertainment model of type ii [expression (20)] with nO = 10, G = 1, or G = 2. As already reported in many articles, increasing κ results in higher proportions of rare alleles in the sample. Plots in Figure 3 also show how ascertainment modifies the distribution of SNP frequencies. Increasing the threshold value G flattens the distribution of frequencies. Both types of ascertainment (i and ii) have similar effects on SNP frequency distributions (results not shown).

Figure 3.

—Probabilities of different configurations of SNP sites, for sample size n = 30; different values of the parameter κ, κ= 0 (○), κ= 1 (*), and κ= 10 (+); under the assumption that data collection was without the ascertainment step [expression (22)] or with ascertainment model type ii [expression (20)] with nO = 10, G = 1, or G = 2.

Likelihood-ratio tests to detect signatures of population growth: An interesting issue is our power to test hypothesis H0 of evolution with constant population size, κ= κ0 = 0, against the alternative hypothesis H1 of population expansion, κ=κ1 > 0, on the basis of SNP data. It is also of interest to determine how this power is affected by the ascertainment step of data collection. From previous computations it follows that SNP data can be seen as samples from multinomial distributions given by expressions (22), (19), or (20). Assuming that the number of SNP sites is always large enough to allow asymptotic approximation (Bickel and Doksum 2001, p. 227) we computed powers of single-value vs. single-value likelihood-ratio tests of statistical null hypothesis H0 (constant population size κ=κ0 = 0) vs. the alternative H1 (population expansion with κ=κ1 > 0). We assumed significance level α= 0.05 and values of κ1, κ1 = 0.1, κ1 = 1, κ1 = 10, κ1 = 100. Table 1, A and B, gives powers of likelihood-ratio tests for sample size n = 50, for different models of ascertainment: no ascertainment [probabilities given by expression (22)] or ascertainment model type ii [expression (19)] with parameters nO and G. Table 1A is for the number of SNP loci M = 30, and Table 1B is for M = 100. From values of powers of tests depicted in Table 1, A and B, one can see that the cases κ0 = 0, κ1 = 0.1 are practically indistinguishable; κ0 = 0, κ1 = 1 may be distinguished only for a large enough number of SNP sites, while κ0 = 0, κ1 = 10, or κ1 = 100 are rather easily distinguishable even for small numbers of SNPs. The ascertainment step in data collection can deteriorate the power to detect signatures of population growth. Increasing the threshold value of G, the aim of which typically is filtering out sequencing errors in the data, also progressively lowers the probability of rejecting the hypothesis H0 of constant population size; i.e., it increases the probability of committing type II error. This results from the flattening effect of increasing G observed in Figure 3.

View this table:
TABLE 1

Powers of likelihood-ratio tests

Data analysis: Data on segregating sites in mitochondrial DNA from Cann et al. (1987): First, we apply our method to the data on segregating sites in mitochondrial DNA from the article by Cann et al. (1987). We fit the exponential scenario (25) to these data by treating each segregating site as an independent SNP. Technically, we estimate the product parameter κ= rNe0 in (25).

Data in Cann et al. (1987) include 195 segregating sites in 148 individuals. Table 2 shows the statistics of segregating sites in these data. Elements in the first column (b) are possible numbers of copies of the rare allele, and elements in the second column (cb) are numbers of segregating sites in the sample that have the number of copies of the rare allele equal to b. Figure 4 shows the plot of log-likelihood function obtained by using expressions (8, 9, 10, 11, 12, 13, 14, 15) and (26). Maximum of the log-likelihood function is attained at κ^=80 . The 95% confidence interval for this estimate, obtained with the use of likelihood-ratio statistics (Bickel and Doksum 2001), is κ ∈ [40, 166].

Segregating sites collected by Cann et al. (1987) are obtained from nonrecombining DNA and the independence assumption is clearly not satisfied. To explore whether a violation of the assumption that SNPs are independent significantly affects the estimate of the parameter κ, we have performed 100 coalescent simulations of genealogies representing ancestries for 148 mtDNA sequences. We added mutations along branches of coalescence trees according to the infinite-sites model with intensity μ. In the simulations we assumed mutational time scale τ= 2μt and exponential change of θ(τ) = 2Ne(τ)μ, θ(τ) =θ0 exp(–ρτ), with parameters θ0 = 400, ρ= 0.2. So, the true value of the product parameter κ was κ= 80. For each of these 100 simulation experiments we treated segregating sites as independent SNPs and we estimated the parameter κ by maximizing likelihood (24). We obtained the mean of estimates equal to 86.8 and standard deviation equal to 29.7. This confirms that our approach, at least for these specific values, will allow us to obtain a reasonable estimate of κ. From these simulations follows the estimate of 95% confidence interval for the parameter κ, when the sample consists of 148 mtDNA sequences and the demography is as shown above. This estimate, [mean – 2 standard deviations, mean + 2 standard deviations], equals κ ∈ [27, 147]. This estimate is quite consistent with the 95% confidence interval obtained from likelihood-ratio statistics, κ ∈ [40, 166]. The shift toward the left of the confidence interval based on simulations results from the asymmetric shape of the distribution of the estimate of κ. By applying a logarithmic transformation to simulation results (estimates of κ) we were able to obtain almost perfect agreement of the two confidence intervals, [40, 166] and [45, 175].

View this table:
TABLE 2

Statistics of segregating sites in mtDNA data

Figure 4.

—Log-likelihood curve for the exponential model of population growth for data on segregating sites in mtDNA from Cann et al. (1987). Each segregating site was treated as a separate SNP. The maximum is attained at Formula .

SNP data from Picoult-Newberg et al. (1999) and Trikka et al. (2002): There are several population studies in the literature where relative frequencies of SNP alleles are shown. We have chosen data from the research by Picoult-Newberg et al. (1999) and data on SNPs in three human genes: BLM, WRN, and RECQL, reported recently by Trikka et al. (2002). In our analysis, we used the data on Caucasians from both sources. The first reason to focus on Caucasians was the possibility of comparing two results, and the second reason was that discovery samples were from Caucasians. Picoult-Newberg et al. (1999, Table 4) have 44 SNP sites in 8 Caucasians (16 chromosomes), while Trikka et al. (2002, Table 2) show allele frequencies of a total number of 31 SNPs in samples of chromosomes of sizes varying from 154 to 158.

When analyzing SNP data we followed remarks given in the source articles (Picoult-Newberget al. 1999; Trikkaet al. 2002) to adjust parameters nA, nO, and G of the model of ascertainment procedure. We modeled the ascertainment procedure for collecting data from Picoult-Newberg et al. (1999) by using expression (17) with nA = 10 and G = 2. A plot of log-likelihood function for the data on Caucasians from the articles by Picoult-Newberg et al. (1999) is shown in Figure 5. It attains a maximum at κ^=3.9 , with the 95% confidence interval, obtained with the use of likelihood-ratio statistics, κ ∈ [0, 105.3]. Log-likelihood function for the data on Caucasians from Trikka et al. (2002) is plotted in Figure 6. Ascertainment was modeled by expression (18) with nO = 10 and G = 1. The maximum-likelihood estimate of the product parameter, from the plot in Figure 6, is κ^=0.78 , with the 95% confidence interval, obtained with the use of likelihood-ratio statistics, κ ∈ [0, 6.1].

Figure 5.

—Log-likelihood curve for the exponential model of population growth for SNP data from Picoult-Newberg et al. (1999). The maximum is attained at Formula .

Sensitivity of estimates to ascertainment model parameters: A question arises: How sensitive are the estimates of parameter κ to changes of the model of the ascertainment? We studied this question by increasing or decreasing the value of the threshold G in expressions (17) and (18). Indexing the estimated parameter with nA, nO, and G, we can denote our estimates from the previous section as κ^[nA=10,G=2]=3.9(Picoult-Newbergetal.1999) (29) and κ^[nO=10,G=1]=0.78(Trikkaetal.2002). (30) Here we compute estimates κ^[nA=10,G=1] , κ^[nA=10,G=3] on the basis of data from Picoult-Newberg et al. (1999) and κ^[nO=10,G=0] , κ^[nO=10,G=2] on the basis of data from Trikka et al. (2002). Analysis of data from Trikka et al. (2002) requires more comment. The model to estimate κ^[nO=10,G=0] assumes that no ascertainment procedure is taken into account. The model to estimate κ^[nO=10,G=2] is inconsistent with complete data of Trikka et al. (2002) in the sense that the data contain one SNP locus with b = 1. To apply the model κ^[nO=10,G=2] we have removed this one locus.

Figure 6.

—Log-likelihood curve for the exponential model of population growth for SNP data from the article by Trikka et al. (2002). The maximum is attained at Formula .

The results of computations show an extreme sensitivity of estimates to the ascertainment model. Notably, κ^[nA=10,G=1]=0,κ^[nA=10,G=3]=(Picoult-Newbergetal.1999) (31) and κ^[nO=10,G=0]=0,κ^[nO=10,G=2]=683(Trikkaetal.2002). (32) In (31), by κ^[nA=10,G=3]= we meant that the likelihood function was increasing for values of κ up to 108.

The fact that the ascertainment model strongly affects estimates of parameters is also confirmed in the previous articles on SNPs. Wakeley et al. (2001) in their Figure 3 show a large bias in SNP frequencies resulting from ascertainment. Similarly, Nielsen (2000)uses a rather oversimplified model, nO = 2, G = 1, for data (Picoult-Newberget al. 1999) and obtains κ^=0 . The need for careful modeling of ascertainment is also stressed by Kuhner et al. (2000).

DISCUSSION

The methods developed in this article allow us to analyze large data sets and carry out computations for different parameter values, which helps us draw more conclusions from data. We have shown examples of applying our methodology to the study of several issues arising in SNP data analysis.

We are particularly interested in the problem of what are reasonable values of the exponential growth product parameter κ= rNe0 obtained on the basis of DNA data. Insight into this problem can be gained by comparing estimates obtained using different approaches.

Our aim when estimating κ from relative frequencies of segregating sites in the article by Cann et al. (1987) was to confirm that different methodologies used for the same data will still lead to comparable results. Therefore, we compared our estimate, obtained by treating nonrecombining segregating sites as SNPs, to those previously obtained on the basis of the same or similar data, but with the use of different methods. Studies that we compared to ours were those by Slatkin and Hudson (1991), Rogers and Harpending (1992), and Polanski et al. (1998), who used pairwise difference statistics, and by Weiss and von Haeseler (1998), who applied the maximum-likelihood approach. Data in these articles originate from different sources, but considering estimations of the authors, reasonable ranges of the product parameter κ, for both the worldwide population and Caucasians, fit into the interval from κ= 50 to κ= 500. Our estimate of κ^=80 is consistent with the above ranges.

Mutation intensity (per site) at autosomal loci is approximately one order of magnitude lower than that in mtDNA (Li 1997). However, the estimate of the product parameter κ= rNe0 is invariant with respect to timescale changes and therefore does not depend on the value of the mutation intensity. We can assume that mutation intensity is used only to scale the time axis. The effective population size for autosomal loci is four times the effective population size for loci at mtDNA. So, the estimate of κ from mtDNA should be one-fourth the estimate of κ from nuclear DNA. Taking into account the large stochastic variation, the estimates of κ coming from SNP data should then be comparable (of the same order of magnitude) to those obtained from mtDNA.

However, our estimates of the parameter κ based on SNP data, κ^=3.9 and κ^=0.78 , are markedly smaller than values coming from mtDNA, which runs counter to the expected tendency. Differences between our estimates and the above ranges can be, probably, attributed to two factors. The first one, mentioned by Wooding and Rogers (2002), is that some fraction of SNPs in the data could be under balancing selection, which would shift their frequencies toward higher values and move the estimate of κ toward lower values. The second factor, which comes from our analysis, is the sensitivity to the parameters of the ascertainment model, shown in (31) and (32). With this high sensitivity, even a small unmodeled factor resulting from eliminating some low-frequency SNPs by assuming that they were sequencing errors can lead to estimates substantially lower than the true value of κ.

Acknowledgments

The authors are grateful to Peter Paule and Markus Schorn for making their program, implementing Zeilberger's algorithm in Mathematica, available to the scientific community. The authors were supported by National Institutes of Health grants GM58545 and CA75432, Polish Scientific Committee (KBN) research projects PBZ/KBN/040/P04/2001 and 4T11F 01824, and NATO collaborative linkage grant LST.CLG.977845.

Footnotes

  • Communicating editor: N. Takahata

  • Received January 29, 2003.
  • Accepted May 30, 2003.

LITERATURE CITED

View Abstract