A two-step method for detecting selection signatures using genetic markers

Summary

A two-step procedure is presented for analysis of θ (F_ST) statistics obtained for a battery of loci, which eventually leads to a clustered structure of values. The first step uses a simple Bayesian model for drawing samples from posterior distributions of θ-parameters, but without constructing Markov chains. This step assigns a weakly informative prior to allelic frequencies and does not make any assumptions about evolutionary models. The second step regards samples from these posterior distributions as ‘data’ and fits a sequence of finite mixture models, with the aim of identifying clusters of θ-statistics. Hopefully, these would reflect different types of processes and would assist in interpreting results. Procedures are illustrated with hypothetical data, and with published allelic frequency data for type II diabetes in three human populations, and for 12 isozyme loci in 12 populations of the argan tree in Morocco.

1. Introduction

The discovery of a massive number of single nucleotide polymorphisms (SNPs) in the genome of several species has enabled exploration of genome-wide signatures of selection via an assessment of variation in marker allele frequencies among populations (e.g. Holsinger & Weir, 2009). Several methods have been proposed for doing this, such as site frequency spectrum, linkage disequilibrium and population differentiation (Sabeti et al., 2006; Akey, 2009). Concerning population differentiation, a parameter θ=F _ST, measuring relatedness between pairs of alleles within a sub-population relative to that in an entire population, has been used for this purpose (Wright, 1951; Cockerham, 1969; Weir & Hill, 2002); Lewontin & Krakauer (1973) and Robertson (1975) discuss related approaches. Equivalently, θ can be interpreted as a measure of dispersion of gene frequencies among groups relative to the variation expected in the population from which such groups derived. For example, Akey et al. (2002) analysed over 26 500 SNPs for which allele frequencies were available in three populations of humans. The θ parameter was estimated for every marker locus and the distribution of estimates over the entire genome, and by chromosome, was examined. By referring these estimates to their empirical genome-wide distribution, 174 candidate genes were identified as possible targets of selection.

Holsinger & Weir (2009) provide an account of the logic of the procedure. Briefly, given a set of loci in a given species, a reasonable assumption is that all share the same demographic history and patterns of migration. If these loci are neutral and have similar mutation rates, members of this set can be conceivably regarded as exchangeable realizations of the same evolutionary process. Loci showing departures from the resulting distribution may serve as flags of genomic regions that have been under the influence of selection. Under the hypothesis of selective neutrality, the distribution (over loci) of estimates of θ is expected to be driven by genetic drift, assumed to affect all loci in a similar fashion. On the other hand, when selection operates on one or several loci (as in a multifactorial model for complex traits), markers that are within genes or in nearby locations will display large or small values of θ, the latter occurring when some sort of balancing selection takes place (Cavalli-Sforza, 1966). This opens an avenue for identification of regions associated with population differentiation, e.g., high versus low producing breeds of dairy cattle. Knowledge of such regions may be useful for enhancing the effectiveness of breeding programs via marker-assisted selection, or for tagging variants associated with disease or quantitative traits. While unusual values of θ may point to genomic locations where selection may have operated, there is arbitrariness with respect to characterizing the type of selection that might have occurred. Typically, loci are classified as either neutral, or subject to balancing selection (low values of θ), or favoured by selection within some specific population or environment (large population differentiation, thus leading to large values of θ). If the values of θ arise from different evolutionary or artificial (such as in plant and animal breeding) processes, one would expect to observe a mixture of distributions leading to clusters representing the different kinds of mechanisms operating. There is no apparent reason why there should be only two or three such clusters; there may be several clusters harbouring loci undergoing different types of selection processes. On the other hand, if θ values vary completely at random due to genetic drift, a single cluster is to be expected.

Statistical issues associated with inferring θ-statistics have been discussed, e.g., by Weir & Cockerham (1984) and Weir & Hill (2002), with emphasis on methods of moments estimation; by Balding (2003) using maximum likelihood for beta-binomial and Dirichlet-multinomial distributions; and by Holsinger (1999), Beaumont & Balding (2004) and Guo et al. (2009) employing Bayesian procedures. None of these treatments have addressed the possible existence of a clustered structure.

The objective of this paper is to present a two-step procedure eventually leading to clusters of θ values. The first step, along the lines of Holsinger (1999), Balding (2003) and Beaumont & Balding (2004), uses a simple Bayesian structure for drawing samples from the posterior distributions of θ-parameters, but without constructing Markov chains. This step assigns a weakly informative prior to allelic frequencies and does not make any assumptions about evolutionary models. The second step regards samples from these posterior distributions as ‘data’ and fits a sequence of finite mixture models, with the aim of identifying clusters of θ-statistics. Hopefully, these would reflect different types of processes and would assist in interpreting results.

The paper is organized as follows. Section 2 reviews basic concepts. In section 3, the first step of the procedure is presented, contrasted with maximum likelihood, and illustrated with a hypothetical dataset and with data on type II diabetes in three populations. Section 4 describes the second step of the procedure, and illustrates it with a dataset containing allelic frequencies for 12 polymorphic isozyme loci in 12 populations of the argan tree (Argania spinosa L. Skeels) of Morocco presented in Petit et al. (1998) and analysed by Holsinger (1999). The paper concludes with a discussion of the proposed methodology.

2. Background

(i) Basic concepts

The stage is set by reviewing essentials of a random effects treatment proposed by Cockerham (1969, 1973). Suppose that genetic markers (e.g. SNPs) are screened in a set of individuals in each of R groups or populations, the latter viewed as drawn at random from some conceptual hyper-population from which such groups derive. Consider a bi-allelic locus (developments carry to multiple alleles as well) and let A _l and a _l be the two alleles at locus l (l=1, 2, …, L); define p _l=Pr (A _l) to be the true, unobserved, frequency of allele A _l in the hyper-population, with 1−p _l=Pr (a _l) being the frequency of a _l. Cockerham (1969) defines a _l as any allele other than A _l and uses an indicator variable x to denote allelic state (‘content’), such that

Here, r=1, 2, …, R denotes group or replicate, i indicates an individual, j is an index for a within-individual deviation and l=1, 2, …, L is an indicator for locus. Even though x _rij,l is a binary variable (so a linear model is questionable), Cockerham (1969) uses the linear decomposition

(1)

where p _l is as before and a _r,l~(0, σ_a,l²), b _ri,l~(0, σ_b,l²) and w _rij,l~(0, σ_w,l²) are mutually uncorrelated zero-mean random deviates, specific to locus l; the σ²′s are variance components. Under the assumption that all alleles at locus l have the same marginal distribution,

and

for l=1, 2, …, L. Decomposition (1) induces the following covariance structure between allelic content variables:

A standard assumption is Cov(a _r, a _r′)=0. The following correlations (all positive) follow.

• • Pairs of alleles drawn at random from different individuals in the same group are correlated as

  (2)

  

  so 0⩽θ_l⩽1 for all l.

• • Pairs of alleles drawn within individuals over all replicates bear a correlation equal to

  

  where F is the total inbreeding coefficient, also known as F _IT (e.g. Weir & Hill, 2002).

• • The correlation between alleles within individuals within the same replicate is

  

  which is the within sub-population inbreeding coefficient f.

It is easy to show that

This expression satisfies

indicating that a reduction in heterozygosity has two sources: one that is due to population sub-division or Wahlund's effect, (1−F _ST,l), and a reduction within subpopulation or group caused by ‘local’ inbreeding, (1−F _IS,l).

Note that parameter F _ST given in (2) can also be written as

In Cockerham (1969), the variance σ_a,l² arises by drawing alleles from a random sample of populations. Under conceptual repeated sampling, this generates a distribution having such variance. However, in many applications, the R groups under study are targeted (as opposed to randomly sampled) populations, e.g., Myles et al. (2007) ; this defines a ‘fixed effects’ model. Now, since σ_a,l² is the between-group variance in allelic content as per model (1), an alternative parametric definition of θ_l in terms of the unknown gene frequencies of the R groups is

(3)

where is the average (over groups) of the frequencies of allele A _l at locus l. Note that _l is taken as an unweighted average; it does not seem sensible to express a parameter in terms of sample size (unless weights assigned to samples reflect true population sizes). Expressing θ_l explicitly in terms of the locus-specific gene frequencies yields

(4)

providing a mapping from the joint space of R allelic frequencies to the single-dimensional space of θ_l, which resides in (0, 1). If allelic frequencies for the different loci are driven from the same stochastic evolutionary process (e.g. as generated by random drift), this defines a distribution of values of θ expected under neutrality assumptions. From a Bayesian perspective, every unknown is a random variable and, since allelic frequencies are unknown, θ as given in (3) will posses a distribution, both a priori and a posteriori. In the first step of the method proposed in this paper, the posterior distribution of θ_l will result from assigning a vague prior to all allelic frequencies, corresponding in some sense to what could be termed as a fixed effects treatment from a frequentist perspective. The second step addresses the question of whether or not all θ_l stem from the same distribution or from different distributions resulting from heterogeneity of the underlying stochastic processes. This makes the treatment proposed here different from those in, e.g., Holsinger (1999) or Balding (2003).

3. Estimation of parameters

(i) Inferring gene frequencies

Gene frequencies can be inferred using a simple Bayesian approach. Suppose that n _r individuals are genotyped in population r, so that the number of alleles screened at locus l is 2n _r=n _r,Al+n _r,al, where n _r,Al and n _r,al are the observed numbers of copies of A _l and a _l, respectively.

A convenient assumption is that of mutual independence between the distributions of alleles at different loci (stronger than that of pairwise linkage equilibrium). Linkage disequilibrium is pervasive but the assumption made above facilitates matters and is widely used, e.g., by Corander et al. (2003). Let p=(p₁, p₂, …, p_R)′ be an RL×1 vector of allelic frequencies for all R groups, where p_r=(p _r,1, p _r,2, …, p _r,L)′ has order L×1. Under the mutual independence assumption, the likelihood conferred by the observed number of copies of alleles to the gene frequencies is

(5)

The maximum likelihood estimator of p _r,l is _r,l=n _r,A _l/2n _r and its empirical variance is . The maximum likelihood estimator is unbiased but unstable, and may take values at the boundaries of the parameter space in small samples.

In a Bayesian treatment, allelic frequencies are assigned a prior distribution that might be homogeneous or heterogeneous over populations, chromosomes or genomic regions (e.g. coding versus non-coding regions). For example, Holsinger (1999, 2006) adopts a prior beta distribution, Beta (p _l|a _l, b _l) (and interprets it as describing variation over populations) with parameters

and

Here θ is common to all loci (i.e. the hypothesis of neutrality) and x _l is the mean allelic frequency at locus l (averaged over populations). Using properties of the beta distribution in the parametric definition of θ leads to

Then, the joint posterior distribution of all unknowns (allelic frequencies, θ and vector x={x _l}) is

Holsinger (1999) took g(θ)=Beta(1, 2) distribution as prior for θ, and assumed that all x _l were identically distributed according to the uniform process g(x _l)=U(0, 1). Given θ and x, the allelic frequencies are conditionally independent with conditional posterior distributions:

where ELSE means all parameters other than p _r,l and the data observed. However, the conditional posterior distributions of θ and x are not recognizable, so an elaborate sampling scheme, e.g., one based on Markov chain Monte Carlo (MCMC) methods, must be tailored. Holsinger (1999) found that inferences were insensitive with respect to the choices of beta and uniform prior distributions for θ and elements of x, respectively. However, it was assumed (as in a neutral model) that all loci share the same θ parameter. This produces a mutual borrowing of information among loci (shrinking p _r,l towards a common value), but the procedures are not explicit with respect to the existence of heterogeneity over loci due to forces such as differential mutation or selective sweeps. As proposed by Beaumont & Balding (2004), one could estimate locus specific θ values and refer these estimates to the posterior distribution of θ under the homogeneity value. In this manner, outliers could be found with respect to the ‘neutral’ distribution, but this would not inform about the structure of any latent heterogeneity.

Here, an alternative approach is used. Jeffreys rule (Bernardo & Smith, 1994; Sorensen & Gianola, 2002) is used to produce a reference prior, which is a distribution assigned to all loci in all populations. This reference prior distribution is minimally informative in a well defined sense (Bernardo and Smith, 1994). Using Bayes theorem, the joint posterior density of all allelic frequencies is now

(6)

Thus, allelic frequencies at different loci are mutually independent, a posteriori, with p _r,l following a beta distribution with parameters and . Possible point estimates of allelic frequencies are the posterior mean

(7)

and the posterior mode

(8)

The variance of the posterior distribution of p _r,l is

(9)

Even though a weakly informative prior is used, differences exist with respect to maximum likelihood. To illustrate this point, consider a hypothetical example with two groups, M and N. Suppose that 100 individuals are genotyped in group M and that the observed number of A _l alleles is 199, i.e. the locus is nearly fixed. The maximum likelihood estimate of p _M,l is 0·995 and its estimated standard error is 4·99×10⁻³; a calculation based on asymptotic normality (without truncation) yields that the probability of obtaining estimates larger than 1 is close to 0·16! Further, the probability of obtaining estimates between 0·9 and 0·995 is close to . On the other hand, the posterior distribution of p _M,l is . The posterior mean and posterior standard deviation are 0·993 (note some shrinkage away from the edge of the parameter space) and 6·06×10⁻³, respectively; the posterior probability of the frequency being larger than 1 is exactly zero, and the probability that p _M,l takes values between 0·9 and 0·995 is about 0·57. Figure 1 displays the posterior distribution of the allelic frequency obtained with Jeffreys prior, overlaid against the normal approximation to the distribution of the maximum likelihood estimates. Clearly, the approach used makes a difference, even in a situation where allelic frequencies are estimated with reasonable precision, as indicated by the small standard error of the maximum likelihood estimate and the small posterior standard deviation in the Bayesian analysis (the coefficient of variation of the posterior distribution is less than 1%).

Fig. 1. Posterior density (thick line) of the allelic frequency p at a locus for which 199 copies have been observed out of 200 alleles counted in hypothetical population M; the posterior distribution is Beta. The thin line is the density of a normal approximation to the sampling distribution of the maximum likelihood estimator.

In the second population, N, 30 individuals are genotyped and 10 alleles are of the type A _l; the maximum likelihood estimate of p _N,l is then , much lower than in M, and its sampling variance is 2·31×10⁻³. The posterior distribution of p _N,l is . In N, the posterior density of p _N,l and the normal approximation to the density of the distribution of the maximum likelihood estimator are very similar (not shown here).

Differences in allelic frequencies between populations M and N at the locus in question may be due to random drift or may suggest a signature of selection.

(ii) Inferring θ by maximum likelihood

A likelihood-based estimate of θ can be obtained by replacing in (3) or (4) the unknown allelic frequencies by their maximum likelihood estimates. For the example of populations M and N above, the estimate is

The sampling variance of the maximum likelihood estimator of θ_l can be approximated using a Taylor series expansion. As shown in Appendix A, the first derivative of θ_l with respect to the allelic frequency at locus l in group r is

for r=1,2, …, R; l=1,2, …, L, where is as before and . Further, let

be an RL×1 vector of first derivatives evaluated at the maximum likelihood estimates of the allelic frequencies. Then, approximately

where is a diagonal matrix containing the estimates of the sampling variances of the maximum likelihood estimates of allelic frequencies p _r,l. For the hypothetical example, . The asymptotic normal approximation to the distribution of the estimates assigns nil probability to ‘estimates’ outside of (0,1); the probability of obtaining estimates of θ between 0·67 and 0·74 for this two-population situation is 0·9996.

(iii) Bayesian inference of θ

Consider now finding the posterior distribution of θ_l as defined in (4) and without making the assumption that the θs are realizations from the same stochastic process, i.e. without borrowing information across loci over and above the shrinkage of allelic frequencies produced by Jeffreys prior. The posterior distribution is analytically difficult to arrive at because θ_l is a non-linear function of gene frequencies in all R groups. However, since it is easy to obtain independent samples from each of the processes, Monte Carlo estimates of features of the posterior distribution of θ_l can be obtained without using MCMC methods at all. Let p _r,l^(s), s=1,2, …, S, be samples from the posterior (beta) distribution of p _r,l, the frequency of allele A _l at locus l. Then, a draw from the posterior distribution of θ_l is given by

(10)

which is a random variable with support in (0,1) (Holsinger, 2006). Then, from S samples, the mean, median, variance, etc., of the posterior distribution of θ_l can be estimated. Each θ_l (l=1, 2, …, L) will have a point estimate and an assessment of uncertainty, e.g., a credibility interval of size 95% given by the 2·5% and 97·5% percentiles of the corresponding posterior distribution estimated either from samples or from the normal theory approximation given in Appendix B.

In the hypothetical populations M and N, the posterior distributions of the frequency of A _l are Beta(199·5,1·5) and Beta(10·5,50·5), respectively. With draws denoted as B ^(s) (.,.), S samples from the posterior distribution of θ_l can be obtained as

To illustrate, 5000 samples were drawn from each of the two beta distributions, to form S=5000 corresponding draws from the posterior distribution of θ_l. The mean and median were 0·6966 and 0·6972, respectively; the standard deviation was 0·070 and the range of values samples spanned from 0·4268 to 0·8883. The posterior density of θ_l and the empirical cumulative distribution function are in Figs 2 and 3, respectively. Values of θ_l appearing with appreciable density range from about 0·5 to 0·9 (Fig. 2), with small posterior probability assigned to values smaller than 0·6 (Fig. 3).

Fig. 2. Posterior density of θ_l for the hypothetical example of populations M and N.

Fig. 3. Empirical cumulative distribution function of θ_l for the hypothetical example of populations M and N.

(iv) A Bayesian ‘null’ distribution for assessing sampling variation uncertainty

It is important to check whether or not posterior estimates of θ_l depart from what would be expected by chance alone. A posterior distribution consistent with expectations under a ‘null’ model is formulated next. The θ_l statistics calculated from the ‘full’ model above can then be referred to this null distribution. Note that the ‘null’ distribution given below describes the uncertainty to be expected from drawing random samples from the same population, but not the variability to be expected due to genetic drift. If estimates of θ_l fall in this null distribution, this would indicate that the study lacks power to answer evolutionary questions in any meaningful manner.

A ‘null’ distribution is arrived at by stating that p _{r , l}=p _l is the same random variable for all R populations. Under this assumption, the posterior distribution of the vector of gene frequencies (now of dimension L×1) under the ‘null’ model is

(11)

Hence, allelic frequencies p _l are mutually independent, a posteriori, with p _l|DATA, Null being a beta distribution with parameters and A draw from the posterior distribution of the F _ST statistic under this model takes the form

(12)

where p _l^(r,s) is a draw from , with R such draws involved in a realization of θ_l^(s), and is the average of the R draws. A set of samples from the posterior distribution of θ_l under the null model is obtained by repeating the sampling process S times. This distribution serves as a reference against which the θ_l statistics calculated from the ‘full’ model can be compared. If the posterior mean of θ_l obtained from the ‘full’ model falls outside of a high density area of the posterior distribution of θ in the null model, then the divergence between populations would be probably due to drift or selection (assuming mutation rates are constant over populations), but not due to chance alone.

For the example of populations M and N, and . Figure 4 depicts the Beta(209·5,51·5) distribution of the allelic frequency under the ‘null’ model. Note that the maximum likelihood estimates of the allelic frequencies in the M and N populations, of 0·995 and respectively, are not assigned any appreciable density under this model. Upon drawing 5000 independent samples from the beta distribution of the allelic frequency under the null model, 5000 draws for θ_l,Null^(s) were obtained by evaluating (12) for each of the samples. Draws ranged from 8·24×10⁻¹³ to 0·0503; the mean (standard deviation) was 0·0038(0·0053) and the posterior median was 0·0017. The posterior density of θ_l,Null was very sharp as shown in Fig. 5. In the full model, the estimated posterior mean (standard deviation) of θ_l was 0·6966, which is unlikely to have been generated under the null distribution. This would make the locus a reasonable candidate for further examination.

Fig. 4. Posterior density of the allelic frequency p under the ‘null’ model for hypothetical populations M and N: 209 copies of A _l are observed out of 260 alleles screened.

Fig. 5. Posterior density of θ_l under the null model for the hypothetical example of populations M and N.

(v) Illustration of sampling variation with candidate genes for type-II diabetes

The Bayesian method was applied to data pertaining to identification of candidate gene variants for type II diabetes in Polynesians (Myles et al., 2007). Prevalence of this disease is high in several Pacific populations, e.g. 40% of adults living in the island of Nauru. DNA samples were obtained from 23 Polynesians, 23 New Guineans and 19 Han Chinese from Beijing. Type II diabetes-associated alleles were from 10 SNP loci having evidence of association. Estimated frequencies and θ_l statistics are shown in page 587 of Myles et al. (2007). To illustrate the Bayesian procedure, data for the KCNJ11 locus was used, and susceptibility allele frequencies (A _l in our notation) were 0·30, 0·25 and 0·34 in the three populations above, respectively. Their figures do not lead to an integer number of alleles, due to rounding error, so the number of observed A _l alleles used here was set to 14 (Polynesians), 12 (New Guineans) and 13 (Han Chinese). Myles et al. (2007) employed an ‘unbiased estimator’ of θ_l for calculating population pairwise differences, and their estimates were 0·003 (New Guinea–China), −0·024 (China–Polynesia) and −0·017 (New Guinea–Polynesia). Note the two negative estimates of a parameter that resides in (0,1); standard errors or significance levels were not provided. Their analysis suggests that this locus is not associated with prevalence of the disease.

The posterior distributions of A _l were Beta_Polynesians(14·5,32·5), Beta_NewGuineans(12·5,34·5) and Beta_HanChinese(13·5,25·5). The number of samples drawn from each of these 3 posterior distributions was S=1000, and 1000 draws from the posterior distribution of θ_KCNJ11 were obtained by evaluation of (10). Values of θ_KCNJ11 ranged from 2·423×10⁻⁵ to 0·1500, with an estimated posterior mean (standard deviation) of 0·019 (0·0193); this estimate is higher than that of Myles et al. (2007). The non-parametric estimate of the posterior density of θ_KCNJ11 is shown in Fig. 6, illustrating that the true value of the F _ST parameter is very likely below 0·10. The posterior inter-correlation structure between allelic frequencies and θ_KCNJ11 in the full model was examined and, as expected, draws from the posterior distributions of allelic frequencies in the three populations were uncorrelated. Samples of θ_KCNJ11 were positively correlated (0·55) with those for allelic frequency in Chinese Han, and the 95% confidence interval for the correlation was 0·51–0·60. However, draws for θ_KCNJ11 were negatively correlated with allele frequencies in Polynesians (−0·07) and New Guineans (−0·39); the confidence intervals for these two correlations were (−0·13,−0·01) and (−0·44,−0·34), respectively.

Fig. 6. Density of the posterior distribution of θ_KCNJ11 obtained from allelic frequencies in Myles et al. (2007).

For the ‘null’ model, the 1000 samples from the posterior distribution of θ_KCNJ11,Null ranged from 3·62×10⁻⁶ to 0·1460, with the posterior mean (standard deviation) estimated at 0·002 (0·002); the posterior median was 0·002 as well. The posterior mean (standard deviation) estimate of θ_KCNJ11 under the ‘full’ model was 0·019, and it did not enter with high density in the ‘null’ model (not shown). Although variation in allelic frequency at locus KCNJ11 among the three populations departs from what would be expected from chance alone (statistical sampling), the observed θ value is very small. This may support the hypothesis that this locus may not be associated with differences in prevalence of type II diabetes, in agreement with Myles et al. (2007). Allelic frequencies were uncorrelated, as it should be, given that the three replicates were drawn from the same Beta(209·5,51·5) distribution. The θ_KCNJ11,Null statistic was uncorrelated with allelic frequencies, and the correlations were −0·08, −0·11 and 0·03 in the three replicates, with all confidence intervals including 0.

4. Clustering of θ-parameters

The second step of the procedure consists of clustering a set of estimates of θ values (in this case, posterior means) from a multi-locus analysis into data driven groups. The expectation is that these clusters might be representative of different processes taking place in the populations such as balancing or directional selection, neutrality or anything else.

The method is illustrated with data from a study of Petit et al. (1998) in which alleles were sampled for 12 isozyme loci of the Argania genus tree in each of 12 areas (populations) of Morocco. The data, given in page 847 of Petit et al. (1998), were modified as shown in Table 1. The modification consisted of treating all loci as bi-allelic by lumping alleles for loci with more than two variants into two classes. The number of individuals sampled per population ranged between 20 and 50, and the number of alleles per locus varied originally between 2 and 5. Note that, at some loci, one of the alleles was fixed in almost all populations. For example, for locus 3, the only population in which segregation was observed was TA.

Table 1. Allelic frequencies at 12 isozyme loci in each of 12 Argan tree populations, adapted from Petit et al. (1998) by making all loci bi-allelic. A1–A12 represent frequencies of the ‘A’ allele at loci 1–12; No. A1–No. A12 are the observed number of copies of the alleles. The number of ‘a’ alleles can be calculated from the number of individuals samples and the number of ‘A’ alleles observed

For each locus, 2000 samples were drawn from the beta posterior distributions of allelic frequencies. For example, the posterior distribution of p _{AB , 1} was Beta(21·5,19·5). From these samples, 2000 draws from the posterior distribution of θ for each locus were formed as in (10). The posterior means were as follows:

so estimates of θ varied over loci from about 0·095 (locus 10) to 0·791 (locus 3); all these estimates did not enter into the corresponding ‘null’ distributions. Box plots of the posterior distributions of the θ parameters are in Fig. 7. Visually, it is tempting to suggest four clusters: the first one would include locus 3, with the posterior mean of θ close to 0·79; the second cluster would include locus 9, with an estimate of θ of 0·59. The third cluster would include loci 5, 7 and 8 with estimates ranging between 0·30 and 0·39, and the fourth cluster would be represented by loci 1, 2, 4, 6, 10, 11 and 12 having the lowest estimates of θ.

Fig. 7. Box plot of the posterior distributions of θ-parameters in 12 isozyme loci of the argan tree in Morocco (data originally from Petit et al. 1998).

The existence of an underlying structure is suggested by the distribution of all 24 000 samples, presented in Fig. 8. In the left panel, a non-parametric density estimate was obtained from these samples treated as if all draws (2000 for each of the 12 loci) had been made from the same process; the densities in the middle and right panels correspond to the logit, i.e., log(θ/(1−θ)), and Gompit, −log(−log(θ)), transforms of the sampled θ values, respectively. The three densities suggest that θ values cluster around 3, perhaps 4, modes.

Fig. 8. Non-parametric density estimates of θ values (based on 2000 sample for each of 12 loci), logit(θ) and Gompit(θ). All samples treated as homogeneous, i.e. as generated from the same stochastic process.

The structure was explored more formally by fitting a sequence of finite mixture models to the means of the posterior distribution of the θ values for each of the 12 loci. These posterior means are independent (under the assumptions made for the allelic frequency models), but not identically distributed, since they are estimated with different precision, due to unequal numbers of individuals sampled and varying allelic frequencies. The distributions of θ values among loci are not normal (the logit and Gompit transforms would be expected to be more nearly so). This should not be an issue because the mixture model was not used for testing hypotheses; its objective, rather, was to explore a clustered structure. Since there are only 12 posterior means, the mixture models must have less than 12 parameters; otherwise, a perfect fit would be obtained. The mixture model fitted to the posterior mean estimates postulated that

where K is the number of components of the mixture (clusters of posterior means of θ values or transforms thereof), π_k is the probability that belongs to cluster k (subject to ), and μ_k and σ_k² are the mean and variance, respectively, of component k. For example, if k=2, there are 5 ‘free’ parameters in the mixture; if k=4, there are 11 such parameters, so it is not sensible to fit a model with more than 4 components. Mixture model parameters were estimated by maximum likelihood via the expectation–maximization algorithm as implemented in the FlexMix package (Leisch, 2004) in the R project (R development core team, 2008). Upon convergence (assuming the stationary point was a global maximum), the conditional probability that (or its transformation) belongs to cluster k is calculated as

The locus was assigned to the cluster with the largest conditional probability. Models with different values of k=1, 2, 3 and 4 were compared using Akaike's information criterion (AIC), that is

where p _K is the number of parameters for a model with K components (McLachlan & Peel, 2000). Models with the smallest AIC values are preferred. It is known that this criterion tends to overstate the number of components due to violation of regularity conditions in mixture models (Celeux & Soromenho, 1996).

Results of the mixture model analysis, by number of components fitted, are shown in Table 2. The AIC favoured a mixture with two clusters when the response was either θ or its Gompit transform, and a single component when the logit transformation was used. Clearly, with data from only 12 loci, the analyses did not have enough power to resolve heterogeneity in a finer manner. This would certainly not be the case with SNP data, where the number of marker loci typically oscillates between a few thousands in some animal species to close to a million in humans. Classification probabilities using K=2 and estimates of cluster mean and standard deviation are shown in Table 3. Irrespective of whether θ values were transformed or not, loci were clustered into two groups, one consisting of loci 3, 5, 7, 8 and 9, possibly reflecting a selection signature, and the other one including the remaining loci, presumably representing neutral loci. The maximum likelihood estimates of the mean and variance of θ values in the cluster with loci 3, 5, 7, 8 and 9 were 0·41±0·21, whereas the corresponding estimates in the other cluster were 0·12±0·03. This assignment into clusters is consistent with the picture emerging from visual consideration of the box plots in Fig. 7.

Table 2. Comparison of mixture models with 2, 3 or 4 components fitted to the 12 posterior means of θ-parameters and their logit or Gompit transforms in the argan tree data of Petit et al. (1998). AIC (models with smallest values are favoured and indicated in boldface)

Table 3. Conditional probabilities of membership to one of two clusters for mixture models fitted to the posterior means of θ for the 12 loci in the argan tree, and their logit, log(θ/1−θ), and Gompit, −log(−log(θ)), transformations (boldfaced probability indicates the cluster with largest probability of membership)

The two-step procedure is simple and does not require the tailoring of problem-specific software. However, it has the drawback of not taking into account the uncertainty associated with the posterior distributions of the θ-parameters, inferred in the first step. In principle, a better approach is to feed the entire set of posterior samples to the clustering procedure, such that not only the location of the posterior distributions of the θs is considered but also their uncertainty as well. Although this is very appealing conceptually, it may create difficulties with the EM algorithm, leading to convergence failure. For instance Qanbari et al. (personal communication) employed the procedure with posterior means (each calculated with 1 million samples from the corresponding posterior distribution) with about 35 000 SNPs in Hereford and Simmental cattle. When posterior means were used as data, the mixture model approach revealed the existence of 4–5 clusters. However, when the 35 million samples were used as data points, the Expectation–Maximization (EM) algorithm, as implemented in FlexMix, failed to converge. An alternative to using the entire collection of samples is to feed a selected set of percentiles of the posterior distribution of each θ_l, so that a proxy for the dispersion of the individual posterior distributions enters into the analysis.

5. Discussion

The use of F-statistics for the study of genetic divergence between population dates back to Wright (1931). Holsinger & Weir (2009) have provided a justification for their usefulness, e.g., in association mapping and in detecting genomic regions affected by evolutionary processes, such as selection. These authors also reviewed different types of statistical methods for inferring F _ST, including Bayesian procedures. Method of moments estimation was prompted by the linear model formalism of Cockerham (1969, 1973), and a review is in Weir & Hill (2002). There has been an increased interest in Bayesian methods, and important contributions in this front have been made by Holsinger (1999, 2006), Beaumont & Balding (2004) and Guo et al. (2009).

In the Bayesian approaches that have been suggested, e.g., Holsinger (1999), the model poses a product binomial (or product multinomial in the case of multiple alleles) likelihood function for allelic frequencies, with conjugate prior distributions, such as beta or Dirichlet processes. Marginalizing over the allelic frequencies yields the beta binomial or Dirichlet-multinomial distributions used by Balding (2003) for likelihood-based inference. Holsinger (1999) matched the mean and variance of, e.g., the beta distribution, to the definition of θ, and obtained a joint posterior distribution which is a function of the unknown allelic frequencies, of θ (assumed exchangeable over all loci) and of the mean allelic frequencies in an undivided population. The implementation, as well as those of Beaumont & Balding (2004) and of Guo et al. (2009) requires MCMC. While the power and flexibility of hierarchical models coupled with MCMC are well known (Sorensen & Gianola, 2002), implementations are not trivial and monitoring of convergence to the equilibrium distribution is a delicate matter (Cowles & Carlin, 1996). The idea in these methods is that, under a neutral model, all θ (over loci) should be realizations of the same stochastic process. Outlying θ values may be suggestive of genomic regions affected by selection. Typically, it is argued that loci are either neutral, subject to balancing selection or to directional selection favouring alleles in specific environments, e.g. Akey et al. (2002). However, the assignment of loci to specific types of processes is often arbitrary.

The present paper follows ideas of Holsinger (1999), but it differs in two important respects. The proposed method has two steps. First, allelic frequencies are assigned a non-informative prior, so that the mutual borrowing of information between loci is limited, leading to less shrinkage of frequencies towards a common value; in maximum likelihood there is no shrinkage at all, an issue criticized by Haldane (1948). Samples of allelic frequencies can be obtained directly (actually, their posterior distributions are tractable, analytically), and these draws are used to form draws from the posterior distribution of locus-specific θ-parameters, using the parametric definition of F _ST as a function of allelic frequencies. The first step was illustrated with hypothetical data and with type II diabetes data in Myles et al. (2007). The step leads to estimates of the posterior distribution of the θs, which can be used to explore underlying structure, presumably caused by different evolutionary forces. In the second step, the structure is explored by using features of the posterior distribution of the θs (posterior means or transformations thereof) as response variables in a mixture model. Data from Petit et al. (1998) on 12 isozyme loci in 12 populations of the argan tree in Morocco were used to illustrate the second step. Here, the posterior means of θ are treated as belonging to a mixture of normal distributions, which is then resolved into data-supported components. Since the final objective is that of clustering loci according to their similarity in θ values, departures from normality are arguably of little consequence. Here, logit and Gompit transformations were examined, and the clustering procedure produced exactly the same results. Using AIC as a gauge for model comparison, it was suggested that the 12 estimates of θ clustered into two groups, one representing putatively neutral loci (provided that this group reflects variation due to drift), and another one possibly corresponding to genomic regions affected by selection. With 12 loci only, it is unreasonable to expect a finer clustering structure. An ongoing study is applying the two-step procedure to large-scale SNP data in an animal population and this will be reported in a future communication.

As mentioned earlier in the paper, the two-step procedure has the disadvantage of not incorporating the uncertainty about the posterior distributions inferred in the first step. Although this can be remedied by using all posterior samples as input into the mixture model analysis, it can create numerical difficulties with the EM algorithm. This is an issue that needs additional research.

The method proposed here extends naturally to multiple alleles. In this case the likelihood is product multinomial, and the beta prior distribution is replaced by a Dirichlet distribution with minimum information content. The posterior distribution of the allelic frequencies is product Dirichlet, which is simple to sample from. Then, samples from the posterior distribution of θ_l would be drawn by evaluation of formulae similar to those in Nei (1973) where θ-parameters are averaged over alleles. For example, one could define

where p _r,l,m is the frequency of allele m at locus l in population r and is the unweighted average over the R populations.

In common with the studies of Holsinger (1999), Beaumont & Balding (2004), Weir et al. (2005) and Guo et al. (2009), the procedure presented here assumes that allelic frequencies are in linkage equilibrium, so that the likelihood of all allelic frequencies is either product binomial or product multinomial. Accommodating linkage disequilibrium, especially with dense batteries of marker loci, represents a formidable task and it is a challenge for future research. For example, Akey et al. (2002) and Weir et al. (2005) reported that θ values of loci in regions of high linkage disequilibrium were similar. Guo et al. (2009) address correlations due to linkage, but not due to linkage disequilibrium, and do so by introducing a spatial structure for loci located in the same chromosome. Specifically, they proposed an autoregressive model in which logit transforms of θ values are correlated according to physical distance. The model is quite involved and requires MCMC computations. However, loci may be in linkage disequilibrium even though not being physically linked (Crow & Kimura, 1970), and such disequilibrium is very common in animal populations (Sandor et al., 2006; de Roos et al., 2008; Lipkin et al., 2009; Qanbari et al., 2009), where finite size and selection under epistasis are factors in building up linkage disequilibrium. The two-step approach considered here could be enhanced by exploring algorithms alternative to EM as well as by consideration of different types of mixtures, e.g., of beta distributions, which are more appropriate for random variables taking values in (0,1).