## Abstract

It is becoming routine to obtain data sets on DNA sequence variation across several thousands of chromosomes, providing unprecedented opportunity to infer the underlying biological and demographic forces. Such data make it vital to study summary statistics that offer enough compression to be tractable, while preserving a great deal of information. One well-studied summary is the site frequency spectrum—the empirical distribution, across segregating sites, of the sample frequency of the derived allele. However, most previous theoretical work has assumed that each site has experienced at most one mutation event in its genealogical history, which becomes less tenable for very large sample sizes. In this work we obtain, in closed form, the predicted frequency spectrum of a site that has experienced at most *two* mutation events, under very general assumptions about the distribution of branch lengths in the underlying coalescent tree. Among other applications, we obtain the frequency spectrum of a triallelic site in a model of historically varying population size. We demonstrate the utility of our formulas in two settings: First, we show that triallelic sites are more sensitive to the parameters of a population that has experienced historical growth, suggesting that they will have use if they can be incorporated into demographic inference. Second, we investigate a recently proposed alternative mechanism of mutation in which the two derived alleles of a triallelic site are created *simultaneously* within a single individual, and we develop a test to determine whether it is responsible for the excess of triallelic sites in the human genome.

BECAUSE of the recent advances in DNA sequencing technologies, it has become feasible to obtain data on sequence variation across tens of thousands of chromosomes (*e.g.*, Coventry *et al.* 2010; Keinan and Clark 2012; Nelson *et al.* 2012; Tennessen *et al.* 2012) and hence to study the impact of variants of very low population frequency. Classical models underlying population genetic studies have typically assumed that each site is affected by at most one mutation event in the genealogical history relating a sample, but for very large samples this assumption is less tenable. One must then account for sites experiencing repeat mutations, which skew the site frequency spectrum and can generate triallelic and even quadra-allelic sites. Triallelic sites are therefore becoming increasingly common, appearing as a few percent of segregating sites in large-scale resequencing studies, particularly as the threshold on masking sites below a given minor allele frequency is being reduced. There are now examples of studies that have found an association between a triallelic single-nucleotide polymorphism (SNP) and a disease phenotype, including coronary heart disease (Crawford *et al.* 2006) and inflammatory bowel disease (discussed in Hüebner *et al.* 2007).

Triallelic sites also have potential use in inference, using frequency spectrum data. The observed frequency spectrum of *diallelic* sites is well recognized as an important summary of genomic data, maintaining a great deal of the information encapsulated by the full data while being relatively simple to interpret. It is therefore well studied: The effects of a host of modeling assumptions on the frequency spectrum have been investigated and many theoretical predictions have been made, typically using either coalescent-based or diffusion-based models. For example, one can obtain analytic results incorporating the effects of a population of varying size (Griffiths and Tavaré 1998; Wooding and Rogers 2002; Polanski and Kimmel 2003; Polanski *et al.* 2003), selection (Griffiths 2003), and population subdivision with instantaneous migration events (Chen 2012). The Poisson random field framework of Sawyer and Hartl (1992) is attractive in this respect because of its amenability to the incorporation of natural selection (Sawyer and Hartl 1992; Bustamante *et al.* 2001). This and other diffusion-based approaches can also be extended to obtain numerical solutions for more complicated underlying population demographic histories, including a single population of variable size (Williamson *et al.* 2005; Evans *et al.* 2007; Boyko *et al.* 2008) or a hierarchy of splitting subpopulations with restricted migration between them (Gutenkunst *et al.* 2009; Lukić *et al.* 2011; Lukić and Hey 2013). Essentially, one writes down the Kolmogorov forward equation for the underlying diffusion approximation and then obtains a numerical solution using finite differences (Williamson *et al.* 2005; Evans *et al.* 2007; Gutenkunst *et al.* 2009) or spectral methods (Lukić *et al.* 2011; Lukić and Hey 2013). Examples of inference using the frequency spectrum such as these are important because they can help us learn about recent human population history, estimate the strength of natural selection, and calibrate our expectations prior to a disease association study. However, none of these approaches make use of the information from triallelic sites since they rely on an infinite-sites assumption in which triallelic sites are never observed (although see Desai and Plotkin 2008; Song and Steinrücken 2012; Steinrücken *et al.* 2013). There have been some extensions to incorporate recurrent mutations into the theory of the frequency spectrum (Sargsyan 2006; Hobolth and Wiuf 2009; Jenkins and Song 2011; Bhaskar *et al.* 2012), but with the exception of Sargsyan (2006) these all assume a simple demography of a stationary, panmictic population of constant size.

In this article, we obtain a closed-form expression for the sample frequency spectrum of a site that has experienced two mutation events, under an extension of the standard coalescent model that allows for very general assumptions about the distribution of times between coalescence events. This allows us to obtain predictions for the shape of the frequency spectrum, allowing for both recurrent mutations and varying historical population size.

To emphasize the usefulness of our results, we consider two applications. First, we investigate the sensitivity of the triallelic frequency spectrum to the assumed demographic history. In particular, our interest is in the following question: How much power to distinguish between different demographic models do we gain by looking at a triallelic, rather than diallelic, site? In a manner quantified further below, we show that although triallelic sites are far less abundant than diallelic sites, they have rather greater value *per site* in capturing the effects of demographic history.

This application relies on a frequency spectrum in which the two mutation events arose independently during the genealogical history of the site. Recently, however, Hodgkinson and Eyre-Walker (2010) noted that there are approximately twice as many triallelic sites in the human genome as would be expected by chance. They explored a number of potential explanations and ultimately favored the idea of a new mutational mechanism: namely, the simultaneous generation of two new alleles due to mutation within a single individual. Although the precise mechanism is unknown, they suggest the instability of base mismatches as a plausible explanation. For example, a mutation of a G = C base pair to an unstable G = A mismatch could give rise to a further mutation to C = A. DNA replication of this mismatch means the ancestral G = C has given rise to both a derived A = T and a derived C = G base pairing (Figure 1). Another possibility is that both strands of the DNA duplex mutate simultaneously due to a chemical or a radiation event. As a second application of our results, we design and implement a frequency spectrum-based test for the hypothesis that a subset of triallelic sites was generated by a *simultaneous* mutation event within a single individual, giving rise to the two derived alleles. The test allows us to account for variable historical population size explicitly, and when we do so we do not find evidence in favor of the existence of such a mechanism (although, as we discuss below, it is likely that this is further confounded by population subdivision in the samples used).

This article is organized as follows. In the following section we first introduce some notation and summarize some previous results. We then obtain closed-form formulas for a triallelic frequency spectrum under a general model of coalescence times distributions. In the sections thereafter we consider two applications. First, we perform an extensive simulation study to discern between the sensitivities of diallelic and triallelic frequency spectra to the underlying demographic model. Second, we obtain the triallelic frequency spectrum under the proposed simultaneous-mutation mechanism of Hodgkinson and Eyre-Walker (2010), and develop a likelihood ratio test to compare it with the null triallelic frequency spectrum under independently-occurring mutations. The test is applied to sequence data taken from the Environmental Genome Project (NIEHS SNPs 2011) and the SeattleSNPs project (SeattleSNPs 2011) to examine whether some fraction of the triallelic sites in these datasets are in fact the product of simultaneous mutations.

## Notation and Previous Results

In this section we introduce our notation and summarize some existing results. Denote by *N*_{0} the diploid effective population size in the present generation and by *u* the probability of a mutation event at a given locus per meiosis. For simplicity we assume throughout that the “locus” is a single site, although we note that the theory extends easily to other loci that may be of interest. Let *θ* = 4*N*_{0}*u* be the population-scaled mutation rate, which we take to be fixed in the usual diffusion limit as *N*_{0} → ∞. In this limit and on a timescale of 2*N*_{0} generations, we denote by *N _{t}* the effective population size at time

*t*back in the past, which we take to be a nonrandom function of time such that

*N*≫ 1 so that a coalescent limit exists for all times (see Slatkin and Hudson 1991 and Griffiths and Tavaré 1994 for details). We assume a general

_{t}*K*-allele mutation model with mutation transition matrix

**P**= (

*P*), so that

_{ij}*P*is the probability forward in time of a mutation taking allele

_{ij}*i*to allele

*j*, given that a type

*i*mutated. It is usual to treat

**P**(and

*K*) as fixed and known. We further denote by

*a*∈ {1, … ,

*K*} the ancestral allele at the site of interest and by

**n**= (

*n*

_{1},

*n*

_{2}, … ,

*n*) the unordered sample configuration taken from that site, with total sample size . A unit

_{K}*K*vector whose

*k*th entry is 1 and all other entries are 0 is denoted by

**e**

*. Finally, let*

_{k}*E*denote the event that there were precisely

_{s}*s*mutation events at the site in the genealogical history relating the sample.

The sample frequency spectrum can be obtained first by finding the probability of the observed sample configuration under the assumptions of an appropriate coalescent model. This may be partitioned according to the number of mutation events in the genealogical history relating the sample. However, for humans the average per-generation mutation rate for SNPs is small; recent studies estimate that on average *u* ≈ 1.2 × 10^{−8} (Campbell *et al.* 2012; Kong *et al.* 2012). Classical population genetics results on the frequency spectrum can be obtained formally by conditioning on precisely one mutation event in the history of the site and then letting *θ* → 0. Denoting the ancestral and derived alleles in a diallelic model respectively by *a* and *b*, it is well known (Watterson 1975; Fu 1995; Griffiths and Tavaré 1998) that for a constant population size, *N _{t}* ≡

*N*

_{0}, and a sample configuration of the form

**n**= (

*n*,

_{a}*n*), (1)since We refer to the quantity

_{b}*φ*(

*i*) as the

*sample*frequency spectrum [as distinguished from the density of the expected number of mutations at each frequency

*x*∈ (0, 1) in a population of genomes comprising many polymorphic sites, which is also referred to as the (site) frequency spectrum]. Throughout this work we obtain the sample frequency spectrum in a finite-alleles model and in the limit as

*θ*→ 0, after conditioning on the required number of mutation events (for triallelic sites, at least two mutation events are of course necessary). In fact, the result (1) is usually obtained by positing a model of

*infinitely many sites*of mutation and then finding the distribution of the number of copies of the mutant allele at any random position at which a mutation occurred. Because we condition on looking at a mutant site, this distribution is equivalent to that of a finite-alleles model at a

*fixed*site and conditioned on one mutation event, with the implicit assumption that

*P*= 0 so that the overall rate of mutation in the two models is the same.

_{aa}There are two extensions to the above result that are relevant to the present work. The first is to general coalescent trees in which the collection of intercoalescence times, **T** = (*T _{n}*,

*T*

_{n}_{−1}, … ,

*T*

_{2}) is not necessarily given by the standard sequence of independent, exponentially distributed random variables. In a standard coalescent model we have that the time

*T*during which there exist

_{k}*k*distinct ancestors to the sample satisfies on the coalescent timescale. However, certain extensions to this model yield a more complicated distribution for

**T**but leave the topological structure of the tree otherwise unchanged. In this setting, Griffiths and Tavaré (1998) have obtained the following result: Under a coalescent model with general intercoalescence times

**T**and conditional on precisely one mutation event at a given site, the sample frequency spectrum is given by (2)where and

*β*=

_{k}*k*. One application of this result is to a coalescent model with a nonconstant population size

*N*. The distribution for

_{t}**T**does not have a simple form, but an expression for (

*T*) is given by Griffiths and Tavaré (1998), and an expression for the marginal density of

_{k}*T*is given by Wooding and Rogers (2002), Polanski

_{k}*et al.*(2003), and Polanski and Kimmel (2003). In the

*Appendix*we provide a new proof of (2), to illustrate our general strategy. For now we merely remark that the topological structure of any polymorphic site having experienced precisely one mutation event in its genealogical history must be of the form shown in Figure 2. Coalescent trees of this form are studied in detail by Wiuf and Donnelly (1999).

A second extension of (1) is to allow for two mutation events at a polymorphic site. In this case we must consider the exact form of the mutation transition matrix **P**. In particular, it may allow for the second mutation to revert a derived allele to its ancestral state (a *back* mutation) or for the second mutation to create a second independent copy of the extant derived allele (a *parallel* mutation). Such mutations do not give rise to triallelic sites, whereas in practice we typically identify sites having experienced two mutations only when three alleles are actually observed. Thus, in extending the definition of the sample frequency spectrum to triallelic sites, we condition on *observing* three alleles, an event we denote *O*_{3}, rather than *E*_{2}. Jenkins and Song (2011) have obtained the following result: Under a standard coalescent model with *N _{t}* ≡

*N*

_{0}the triallelic sample frequency spectrum is given by (3)whereandIn the above expression

*a*,

*b*, and

*c*are distinct alleles with

*n*+

_{a}*n*+

_{b}*n*=

_{c}*n*;

*a*is the ancestral allele and

*b*and

*c*are derived alleles. We set

*H*

_{0}:= 0 by convention. If the diagonal of

**P**is zero, then the sums involving

**P**in

*C*respectively simplify to [1 − (

**P**

^{2})

*] and [1 − (*

_{aa}**PP**

*)*

^{T}*]. (That this simplification requires the diagonal of*

_{aa}**P**to be 0 was inadvertently omitted from Jenkins and Song 2011, Corollary 6.1.)

## General Triallelic Frequency Spectrum

In this section we obtain a closed-form expression for the sample frequency spectrum of a triallelic site under a general coalescent model with variable population size, *N _{t}*. This generalizes (2) to the case of two mutation events at a single site and generalizes (3) to the case of a variable population size. Our arguments and notation are similar to those in Jenkins and Song (2011), and a brief proof is deferred to an Appendix. In that article, a key observation is that the event

*E*

_{2}can be partitioned as follows:

*E*_{2}: The two mutation events are genealogically nested.*E*_{2}: The two mutation events are genealogically nonnested and at least one of them does not reside on the basal (adjacent to the root) branches of the tree.*E*_{2ℬ}: The two mutation events reside on the two different basal (adjacent to the root) branches of the tree.*E*_{2S}: The two mutation events reside on the same branch of the tree.

Only the first two cases can lead to a triallelic site, and so we do not consider the last two any further. It is straightforward to obtain analogous generalizations for the last two cases, although we omit them. The events *E*_{2} and *E*_{2} are illustrated in Figure 3. In these examples, the older of the two mutation events gives rise to the allele *b* and the younger gives rise to the allele *c*, and the subsets of *E*_{2} and *E*_{2} satisfying these constraints are denoted by and respectively. To find the sample frequency spectrum we first consider the joint probability of observing our triallelic sample with each of these events.

**Lemma 1.** *Let n _{a}*

**e**

*+*

_{a}*n*

_{b}**e**

*+*

_{b}*n*

_{c}**e**

*(3)*

_{c}denote a triallelic sample as in*. Under a coalescent model with time-dependent population size N*,

_{t}*the joint probability of such a sample together with the way the two mutations are placed on the coalescent tree satisfies*(4) (5)

*as θ*→ 0

*. Furthermore*, (6) (7)

*The coefficients in the above expressions are*

*where δ*

_{j}_{,}

*.*

_{k}denotes the Kronecker delta*Proof*. See the *Appendix*. □

From Lemma 1, we can obtain our main result in a straightforward manner.

**Theorem 1.** *Let n _{a}*

**e**

_{a}+ n_{b}**e**

_{b}+ n_{c}**e**

*:*

_{c}denote a triallelic sample as above*a*,

*b*,

*and c are distinct alleles with n*;

_{a}+ n_{b}+ n_{c}= n*a is the ancestral allele and b and c are derived alleles. Under a coalescent model with time-dependent population size N*,

_{t}and in which mutation events occur independently in the tree*the sample frequency spectrum is*(8)

*where*

*and*{⋅}

*denotes the indicator function*.

*Proof*. As in Jenkins and Song (2011, Theorem 6.2), this follows fromNow substitute for each term on the right-hand side using Lemma 1 and let *θ* → 0.□

Thus, while the frequency spectrum for a site experiencing one mutation event depends only on the first moments of the elements of **T** [see (2)], the frequency spectrum for a site experiencing two mutation events depends only on the second moments of the elements of **T** (Theorem 1). These moments are considered in further detail by Polanski *et al.* (2003) and Zˇivković and Wiehe (2008). Under a suitable choice of historical population size function, *N _{t}*, the frequency spectrum given by (8) serves as our null model for triallelic sites. As a check on (8), we can fix the population size,

*N*≡

_{t}*N*

_{0}, so thatInserting this expression into (8) leads to (3), after extensive simplification.

While this article was under review we learned of related work by Sargsyan (2006), who also obtains an expression for the frequency spectrum of a site experiencing two mutation events under an arbitrary distribution on **T** (Sargsyan 2006, Lemma 34). Our work strengthens his result, which relies on higher-order and exponential moments of the elements of **T**. Our work also allows for a more general model of mutation and disentangles the relative contributions of nested and nonnested mutations.

## Application I: Sensitivity to Demography

In this section, we compare the frequency spectrum of a diallelic site with that of a triallelic site. Given a sample taken from a population whose recent history is described by a demographic model , we can measure the information that is lost if one erroneously applies the frequency spectrum according to another model . To quantify this difference in information, we employ the Kullback–Leibler (KL) divergence, a measure defined for this task (Kullback and Leibler 1951; Burnham and Anderson 2002). We define the KL divergence from to bywhere is the appropriate sampling distribution under model and denotes expectation with respect to random samples **n** drawn under model . Thus, KL divergence is the expected likelihood ratio when testing an alternative against a null and the alternative is true. Although KL divergence properly refers to distributions under these models rather than to the models themselves, when we refer to the divergence between two models, it should be clear from the context that we are referring to either their diallelic or their triallelic sample frequency spectrum.

We focus on the divergence from one model of population growth to another. The KL divergence (amount of information loss) can thus be compared for samples from a diallelic site *vs.* samples from a triallelic site. A larger value of KL divergence for triallelic sites would suggest that such sites are potentially very informative for demographic inference. Throughout this section, a symmetric mutation matrix is used in the frequency spectra calculations (*i.e.*, it is assumed that all transitions between alleles are equally likely). To illustrate how divergences vary at different scales, we focus our analysis on frequency spectra for samples of 10 and 100 individuals [a typical magnitude of sample sizes in demographic inference studies (Williamson *et al.* 2005; Gutenkunst *et al.* 2009)].

### The effect of sample size in simulations

To compute frequency spectra under general models of historical population size we need first- and second-order moments of **T** (*cf*. Equation 2 and Theorem 1). We precompute these by simulating coalescent trees, using ms (Hudson 2002). To investigate the effect of sample size in simulations and of the differing dimensions of the two frequency spectra, we first consider the case of a constant population size (*i.e.*, *N _{t}* ≡

*N*

_{0}), for which the expected frequency spectra are known in closed form (see Equations 1 and 3). Specifically, we compute the KL divergence from the expected frequency spectrum computed exactly to the expected frequency spectrum obtained by simulation of

*N*

_{trees}to approximate the first- and second-order moments of

**T**. The results are shown in Figure 4. Although it might be considered unfair to compare KL divergences in diallelic frequency spectra (one-dimensional distributions) with those between triallelic spectra (two-dimensional distributions), Figure 4 clearly shows that these divergences exhibit extremely similar behavior as we increase

*N*

_{trees}. The degree to which the true spectrum is approximated by its Monte Carlo counterpart is almost exactly the same in the di- and triallelic cases, regardless of the choice of

*N*

_{trees}. We also note that for both di- and triallelic spectra, the KL divergences have similar magnitude between sample sizes 10 and 100 for all choices of

*N*

_{trees}. Thus, KL divergence appears to be a good measure of the difference between two frequency spectra that is relatively invariant to the differences in dimensionality and sample size in our study. Because Figure 4 illustrates that using

*N*

_{trees}= 10

^{6}results in negligibly small divergences (on the order of 10

^{−8}) from the true (closed-form) frequency spectra, we fix

*N*

_{trees}= 10

^{6}in the remainder of this section.

### Exponential growth

Next, we examine how sampling from a population with historical exponential growth affects the resulting frequency spectrum. As specified in Figure 5, we investigate seven models of exponential growth, * _{i}*,

*i*= 1, … , 7 (with

_{0}representing a population of fixed size). To compute their respective moments of intercoalescence times, we first simulate 10

^{6}trees from populations according to each model. In Figure 6, A and B, we compute the KL divergence D(

*‖*

_{i}_{0}) of the sample frequency spectrum under

_{0}from the sample frequency spectrum under

*, for*

_{i}*i*= 1, … , 7. To investigate the potential benefit of triallelic spectra in fine-tuning between two population growth models with different degrees of exponential growth, we also compute KL divergences

*D*(

*‖*

_{i}

_{i}_{−1}) of sample frequency spectra under growth model

_{i}_{−1}from growth model

*, for*

_{i}*i*= 1, … , 7 (Figure 6, C and D).

Figure 6 demonstrates clear superiority of using triallelic spectra to distinguish between demographic models with varying degrees of exponential population growth. The mean KL divergence from exponential growth model *i* to *i* − 1 is increased by 87% when triallelic spectra are used in place of diallelic spectra for samples of size 10 (and the divergence is increased by 99% for sample size 100). This indicates that triallelic sites contain information that may significantly increase our ability to discern between competing exponential growth models with similar parameters.

### Instantaneous growth

We next investigate in further detail the effect of sample size on KL divergence from a growth model to a model of fixed population size. For the growth model we assumed a function of historical human population growth as inferred by Williamson *et al.* (2005), who assumed an instantaneous expansion of the population from an ancestral size *νN*_{0} to a modern size *N*_{0} a time *τ* ago. Using data from the Environmental Genome Project and working within the framework of the Poisson random field model (Sawyer and Hartl 1992), Williamson *et al.* (2005) inferred maximum-likelihood estimates (MLEs) of and the latter in units of 2*N*_{0} generations. (The authors estimated *N*_{0} ≈ 51,340 directly by comparing polymorphism and divergence data, in which case the latter estimate is calibrated as generations or, further assuming a 20-year generation time, years.) We denote this model by _{W}. Thus, given samples that actually stem from a population such as the one described by Williamson *et al.* (2005), the KL divergence *D*(_{W} ‖ _{0}) quantifies the ability to distinguish that these samples do not come from a fixed-size population. This analysis therefore studies the effect of sample size on the ability to perform inference of population growth parameters under realistic settings in a problem of great interest, for both di- and triallelic sites.

Figure 7 illustrates that the KL divergence using triallelic spectra for the two models is much greater, for any sample size, than the divergence using the corresponding diallelic spectra. Furthermore, we see that the increase in KL divergence that results from the presence of the third allele at a triallelic site also grows with increasing sample size, at least up to sample sizes of 25, before leveling off beyond 25 and providing a consistent ∼97% increase in KL divergence.

We further address the effect of the parameters of a model of instantaneous population size change, as follows. Adopting _{W} as a reference point, we examine variations in the two parameters *ν* and *τ*. First, the amount of instantaneous growth is varied while keeping the time of the size change fixed to same value as _{W} (see models 8–20 in Table 1), and subsequently, different times of size-change occurrence are examined while the amount of instantaneous growth is fixed so that the prechange population size is 15% of the postchange size (models 21–29 in Table 1).

Repeating the steps of our analysis of the exponential growth models, we computed KL divergences *D*(* _{i}* ‖

_{0}) for

*i*= 8, … , 29 (Figure 8, A and B, and Figure 9, A and B), and to investigate the potential benefit of triallelic spectra in the more subtle problem of distinguishing between two instantaneous population growth models with different degrees of growth, we computed

*D*(

*‖*

_{i}

_{i}_{−1}) for each

*i*= 9, … , 20 and

*i*= 22, … , 29 (Figure 8, C and D, and Figure 9, C and D).

From Figure 8 and Figure 9, we again find that KL divergence from the spectrum under one growth model to another is larger when we use a triallelic, rather than diallelic, sample frequency spectrum, although this advantage fades as the time of the sudden size change is moved extremely far into the past while the amount of instantaneous growth is kept constant (as in models 21 and 22 in Figure 9). The advantage of the triallelic sample frequency spectrum is even more pronounced for sample size 100 than for sample size 10. The mean KL divergence *D*(* _{i}* ‖

_{i}_{−1}) (over

*i*= 8, … , 20; see Figure 8, C and D) is increased by 79% when triallelic spectra are used in place of diallelic spectra for samples of size 10 (and this divergence is increased by 97% for sample size 100). For models 21–29 (Figure 9, C and D, the mean KL divergence

*D*(

*‖*

_{i}

_{i}_{−1}) is increased by 56% through the inclusion of the third allele for samples of size 10 (and the divergence is increased by 83% for sample size 100). Thus, the inclusion of the third allele in triallelic spectra contains information that may considerably increase our power to discern between competing instantaneous growth models with similar parameters.

## Application II: Simultaneous Mutation Model

### Theory

As discussed above, Hodgkinson and Eyre-Walker (2010) propose that there exists another mechanism of mutation responsible for the observed excess of triallelic sites in samples of human genomes: the simultaneous generation of two new alleles within a single individual. It is estimated that this mechanism is responsible for the generation of ∼3% of all human SNPs. In support of this hypothesis they developed a phylogenetic statistic to test whether the two minor alleles of a triallelic site are closer to each other on a reconstructed phylogenetic tree than would be expected by chance. Using this test, they find significant evidence at the 5% level for proximity of the minor alleles when probabilities are combined across all triallelic sites in their data, although the null hypothesis of two independent mutation events is rejected at a rate close to the nominal 5% when each of 113 triallelic sites is tested independently. There are, however, some limitations to this test. First, as Hodgkinson and Eyre-Walker (2010) observe, phylogenies are reconstructed using local haplotype information but ignoring the confounding effects of recombination. Second, it uses the mean branch length between leaf nodes subtended by minor alleles as an indirect measure of the branch length between the minor alleles themselves. Finally, it uses the *proximity* of two mutation events on the phylogeny as evidence for what is in fact a stricter hypothesis of *colocation*. In this section we use our results on the triallelic frequency spectrum to develop a complementary test that does not suffer from these issues. Our approach is to find the sample frequency spectrum as predicted by the simultaneous mutation mechanism and then to compare it with the results of Theorem 1 via a likelihood-ratio test. It will be clear that the two spectra give very different predictions, particularly with regard to the expected number of singleton alleles in a sample.

The key observation that enables us to obtain the frequency spectrum under this model is as follows. Suppose there exists a mechanism whereby two new alleles are produced within a single individual, such as that described in the Introduction: A single DNA duplex within a diploid cell experiences subsequent mutations of both of the nucleotides within a base pair. The duplex then undergoes replication so that two new alleles are produced (Figure 1). Now, the two alleles are observed in a sample taken from the population in the present day. The individual responsible for the creation of these alleles must have been an ancestor of individuals in the sample carrying *either* of the derived alleles. Moreover, this ancestor was the *most recent common ancestor* of any pair of individuals carrying the two distinct derived alleles. The genealogy relating the sample at this site must be of the form illustrated in Figure 10; in particular, the simultaneous mutation event coincides with the coalescence node uniting the two clades defined by individuals carrying the two derived alleles. We can therefore condition on this coincidence event in deriving the sample frequency spectrum under this model, by choosing uniformly among the *n* − 1 coalescence nodes. As noted by Hodgkinson and Eyre-Walker (2010), the probability that both products of a single human meiosis leave descendants in the following generation is negligible, so the posited simultaneous mutation event is presumed to occur during the mitotic phase of germ-line development.

In its full generality, the mutations under this model are parameterized by a transition matrix **Q** = (*Q _{i}*

_{,{}

_{j,k}_{}}) whose (

*i*,{

*j*,

*k*})th entry specifies the probability that a simultaneous mutation affecting allele

*i*gives rise to alleles

*j*and

*k*. For notational simplicity we assume

*Q*

_{i}_{,{}

_{j}_{,}

_{k}_{}}= 0 if

*i*,

*j*,

*k*are not all distinct—it is straightforward to make the appropriate modifications to relax such an assumption. In this setting we have the following theorem.

**Theorem 2.** *Let n _{a}*

**e**

_{a}+ n_{b}**e**

_{b}+ n_{c}**e**

*,*

_{c}denote a triallelic sample*and let*

*denote the event that there were s instances of the mechanism of simultaneous mutation in the genealogical history relating the sample. Then the sample frequency spectrum is*(9)

*Proof*. See the

*Appendix*.□

We remark that in the above we conditioned on ; equivalently one could introduce a rate parameter for the occurrences of simultaneous mutations and then let it go to zero after conditioning on *O*_{3} only, in which case none of , , … contributes to the frequency spectrum.

Importantly, the frequency spectrum in (9) depends on the distribution of topologies of coalescent trees but not on the distribution of **T**. Thus, Theorem 2 continues to apply when we allow a general distribution of intercoalescence times as in the null model, and in particular this includes a model of variable population size, *N _{t}*. To summarize: Under a model in which the population size

*N*is allowed to vary in time, the sample frequency spectrum of a triallelic site is given by (8) when the two mutation events occur independently and by (9) when they occur simultaneously within a single individual. An example of the two spectra is shown in Figure 11. Clearly, the largest difference occurs in the frequency class of double singletons, (

_{t}*n*,

_{a}*n*,

_{b}*n*) = (

_{c}*n*− 2, 1, 1), which contributes 0.37 of the total probability mass under the simultaneous mutation mechanism compared with 0.09 when the two mutations occur independently along the tree. Other nearby configurations in which both derived alleles are at very low frequency are also overrepresented according to the simultaneous mutation mechanism by comparison with independent mutations, while configurations in which one or both derived alleles are at moderate frequency are slightly underrepresented. This underrepresentation is greatest for frequencies of the form (

*n*, 1,

_{a}*n*) and (

_{c}*n*,

_{a}*n*, 1),

_{b}*i.e.*, along the axes in Figure 11.

### Likelihood-ratio test of independent *vs.* simultaneous mutation

Our goal now is to test for a relative excess of triallelic sites that conform to the frequency spectrum of the simultaneous mutation mechanism. We take as our null hypothesis that each triallelic site was generated by two independent mutation events, so that the frequency spectrum is given by *φ*_{0} (Equation 8). An appropriate alternative is that some fraction, *λ* > 0, of triallelic sites arose as the result of a simultaneous mutation event. Under this model, the sample frequency spectrum is given by the mixture (10)Suppose we observe *M* triallelic sites, and the *i*th site has configuration If each pair of sites is sufficiently far apart that their genealogical histories are independent, then a likelihood-ratio statistic for these data is (11)where is a MLE for *λ*, with We reject the null hypothesis at level *α* if Λ lies within the 100(1 − *α*)th percentile tail of its null distribution. Unfortunately, the null distribution of −2 ln Λ does not tend toward the usual -distribution due to the possibility that the mixture parameter *λ* lies on the boundary of its permissible set (Self and Liang 1987). We thus employ bootstrap estimation to determine the null distribution of the test statistic by simulation (further described in supporting information, File S1).

Our aim is to apply this test to empirical SNP data, so we take *K* = 4 with alleles {A, C, G, T}. However, to apply the test we must also specify both transition matrices **P** and **Q**. While **P** can be estimated by, for example, using the empirical frequencies of each type of mutation event inferred from diallelic SNPs and their corresponding outgroup alleles (Chan *et al.* 2012), there is no guidance on how to choose **Q**. In fact, we do not expect the test to be greatly influenced by the *identities* of the three alleles at a site compared to the information contained in the sample counts themselves. Therefore we choose to eliminate the appearance of the entries of **Q** in Λ by *conditioning* on the identities of the observed alleles at each site. Formally, we replace *O*_{3} in the definitions above with , the event that three alleles are observed and these alleles are *a*, *b*, and *c*, with *a* ancestral. This leads to the slightly modified test statistic , obtained by replacing the appearances of *φ*_{S} and *φ*_{0} in each of (10) and (11) withBy repeating the reasoning that led to expressions for *φ*_{S} and *φ*_{0}, we obtain from *φ*_{S} by dividing by *Q _{a}*

_{,{}

_{b}_{,}

_{c}_{}}in (9) and from

*φ*

_{0}by replacing

*κ*

_{j}_{,}

*in (8) withThe resulting test statistic is independent of*

_{k}**Q**and is more robust than Λ to the choice of

**P**.

### Data

To apply our test to SNP data we first need to specify **P**, and to do this we followed the procedure described in Chan *et al.* (2012). Their method requires empirical counts of each type of diallelic SNP with the ancestral allele at each SNP specified, as well as the overall abundance of each nucleotide in the genome. To obtain empirical diallelic SNP data we used the Genome Variation Server (v6.01) (GVS 2011). We obtained each diallelic SNP from the GVS database, discarding those for which the orthologous chimpanzee allele was unavailable or did not match any of the human alleles. To restrict our attention to neutral mutation events occurring at a typical genomic rate, we further discarded SNPs at which one of the alleles would produce a CpG dinucleotide; we discarded SNPs residing in coding regions; and to keep our estimate of **P** independent of the test data set we discarded sites at which more than two alleles were observed. Assuming that the chimp carries the ancestral allele at each site polymorphic in humans and that each remaining SNP represents a single mutation event from the ancestral to the derived allele, we were left with the following counts of each type of mutation event:These counts can be converted to rates of mutation by comparison with the overall genomic abundance of each type of nucleotide, which were obtained from the University of California, Santa Cruz Genome Browser (hg19) (Kent *et al.* 2002), after excluding both CpG dinucleotides and the Y chromosome and weighting the counts for the X chromosome by 3/4. This left the following counts of each type of nucleotide in the human genome: *N _{A}* = 835,878,173 (30.1%),

*N*= 836,874,687 (30.0%),

_{T}*N*= 552,795,868 (19.9%), and

_{C}*N*= 553,090,147 (19.9%). Using these values together with the mutation counts given above, we obtained the following empirical

_{G}**P**matrix, using the method described in Chan

*et al.*(2012):

We reanalyzed the triallelic data set of Hodgkinson and Eyre-Walker (2010, their Table S1), which comprised 113 triallelic sites. These were in turn obtained from 896 nuclear genes sequenced as part of the Environmental Genome Project (NIEHS SNPs 2011) and the SeattleSNPs project (SeattleSNPs 2011), which provide high-quality resequencing data, avoiding problems such as ascertainment bias. Only sites of high quality (*Q* > 25), outside CpG dinucleotides, and outside coding regions were included in the data. Orthologous chimpanzee alleles corresponding to each site in the data were kindly provided to us by Alan Hodgkinson, CHU Sainte-Justine, Montréal, Canada; sites for which the chimp allele was unavailable were excluded from our analysis, leaving *M* = 96 triallelic sites. Since these sites originate from different experiments using different population panels, the sample size varied across sites. The minimum, mean, and maximum sample sizes across the 96 triallelic sites were 71, 160, and 190, respectively.

To compute Λ, the required formulas are given in terms of the joint moments [*T _{j}T_{k}*], and as before we precompute these numerically. This precomputation step can be reused for different choices of

*λ*and across segregating sites with the same sample size, and so it does not add to the computational burden significantly.

The most striking feature of the historical human population size is its recent rapid growth (Keinan and Clark 2012). To examine this, we used the model _{W} of Williamson *et al.* (2005) as described above. Computing the MLE under this demographic model yielded conflicting with the conclusions of Hodgkinson and Eyre-Walker (2010). This illustrates the importance of accounting for demographic changes when using frequency spectrum data; had we assumed a population of constant size we would find (*P* < 0.001; see File S1, in which we also examine the robustness of this result to assumptions about **P** and to potential sequencing errors), in better agreement with Hodgkinson and Eyre-Walker’s (2010) estimate of *λ* ≈ 0.5. The decrease in the value of our estimated mixture parameter can be explained by the fact that singleton sites—and, similarly, doubly singleton triallelic sites—are relatively more probable under a null model with population growth than one without. Thus, the abundance of doubly singleton sites that previously gave rise to an extreme likelihood-ratio statistic is now explicable under the null model.

A further potential complication of the data used here is population subdivision with migration between subpopulations. Analytic results for the frequency spectrum under complex demography are unavailable, even ignoring the issue of recurrent mutations (but for recent progress on this problem see Chen 2012), and so in File S1 we further investigate this issue by a simulation study. We find although the lack of availability of subpopulation labels with the data renders the result nonsignificant (*P* ≈ 0.42; File S1).

## Discussion

In this article we have obtained closed-form expressions for the frequency spectrum of a site experiencing two mutation events, of which triallelic sites are an important example, and have allowed for the possibility of a varying historical population size; the results generalize those of Griffiths and Tavaré (1998) and Jenkins and Song (2011). We applied our formulas to the question of the ability of the frequency spectrum to discern between closely related models of population growth and to the question of the mechanism of mutation that gives rise to triallelic sites.

Demographic inference from SNP data has thus far relied solely on diallelic sites. As sample sizes in sequencing studies grow with the falling cost of the technology, it is likely that an ever-increasing fraction of segregating sites found will be triallelic. In this article, we have illustrated that the triallelic sampling frequency spectrum is more sensitive to historical population size changes under both exponential and instantaneous growth models, and thus it seems likely that improved estimates of population growth parameters may be obtained by incorporating this growing number of triallelic sites in demographic inference analyses. Furthermore, the increased sensitivity of the triallelic spectrum over the diallelic spectrum to discern between different growth models becomes more exaggerated as sample sizes are increased. While triallelic sites remain relatively rare compared to diallelic ones, the above analysis suggests that they are more valuable per site in distinguishing between a variety of demographic models. We hope that the ability of triallelic sites to fine-tune between competing growth models will be especially useful when looking at recent *superexponential* growth (Keinan and Clark 2012), although convenient software for this type of growth is not yet available.

We also found the frequency spectrum under a model in which the two derived alleles of a triallelic site may be generated simultaneously (Hodgkinson and Eyre-Walker 2010), and we have developed a likelihood-ratio test for the existence of such a mechanism. This approach is parameterized by the mixture parameter *λ* that represents the fraction of triallelic sites having arisen as a result of the simultaneous mutation mechanism. Assuming a simple randomly mating population of constant size, we find a MLE of which is significantly nonzero. We show that another explanation for the excess of doubly singleton triallelic sites is rapid recent population growth, but when we posit a realistic demographic model that includes population subdivision and migration as well as recent growth, we find only a minor adjustment to supporting the idea that at least some triallelic sites were generated as the result of a simultaneous mutation event. This latter estimate is not, however, significantly different from 0, a state of affairs we can at least partly attribute to a considerable loss of power—the individuals sampled at a substantial fraction of triallelic sites in our data set are lacking subpopulation labels. [In particular, none of the triallelic sites of the form (*n*^{(}^{i}^{)} − 2, 1, 1) remained triallelic after removing from the sample individuals of unknown origin.]

Because of the lack of power associated with the available data, we treat these results with caution and do not rule out the possibility that at least some sites were generated by a mechanism of simultaneous mutation. Indeed, most of the information about a mechanism generating an excessive number of triallelic sites is contained in the absolute *number* of triallelic sites observed in the human genome. We did not use this quantity; instead we *conditioned* on the observed number *M* of triallelic sites and addressed a slightly different question: Given the excessive number of triallelic sites in the genome, is the frequency spectrum of some fraction of these sites consistent with the two derived alleles being generated simultaneously within a single individual? Even if this proposition is rejected, the question of why there is such an excess of triallelic sites remains. We are hopeful that the coming flood of (subpopulation-labeled) genomic data will enable us to address these questions with much improved power.

One lesson of our work is that rare alleles are as vital when looking at triallelic sites as when looking at diallelic ones; singleton alleles at *diallelic* sites are already the focus of other coalescent-based tests of neutrality (Fu and Li 1993; Achaz 2009). While earlier genotyping projects often chose to exclude sites below a given minor allele frequency, typically 5%, more recent trends—improved sequencing technologies, larger sample sizes, and an interest in rare variants in their own right (Cirulli and Goldstein 2010; Coventry *et al.* 2010; Keinan and Clark 2012; Nelson *et al.* 2012; Tennessen *et al.* 2012)—should serve to make the appropriate data more readily available. As this trend continues, we expect our results to find further use as recurrent mutations manifest themselves more and more commonly.

## Acknowledgments

We thank Christoph Theunert for helpful discussions and Alan Hodgkinson for providing ancestral allele information to accompany the triallelic site data. We also thank John Wakeley for helpful suggestions for improving the exposition of this article. This research is supported in part by National Institutes of Health grant R01-GM094402, an Alfred P. Sloan Research Fellowship, and a Packard Fellowship for Science and Engineering.

## Appendix

*Proof of* (2). The arguments given here mimic in part those found in Jenkins and Song (2011), and we refer the reader to that work for further details. In that article we work with *unordered* configurations, **n**; here, the argument is easier to illustrate using *ordered* configurations. We denote a random vector of *n* alleles consistent with the unordered configuration **n** by **v _{n}**; by sampling exchangeability there are equiprobable such vectors. First, denote the event that a single mutation occurred in the genealogical history relating the sample and that it gave rise to a derived allele

*b*by . Now, inspection of Figure 2 tells us that any particular coalescent history with leaf configuration

**v**and consistent with must exhibit the following sequence of events going back in time, for some

_{n}*l*∈ {1, 2, … ,

_{a}*n*}: a collection of

_{a}*n*−

_{a}*l*coalescence events of type

_{a}*a*lineages and a collection of

*n*− 1 coalescence events of type

_{b}*b*lineages, with the two collections in some interspersed ordering; a mutation event taking the single remaining type

*b*lineage to a type

*a*lineage; and coalescence of the remaining

*l*+ 1 type

_{a}*a*lineages to a most recent common ancestor of the sample.

The quantity *l _{a}* represents the number of extant lineages ancestral to samples of type

*a*at the time of the single mutation event. Suppose the first coalescence event is between two lineages whose allele at the leaves of the tree is

*b*; such a coalescence occurs (among all possible coalescence events) with probabilityContinuing in this vein back to the most recent common ancestor, we find that the probability of a given sequence of

*n*− 1 alleles corresponding to each coalescence event, which is consistent with

**v**and with satisfiesFor example, in Figure 2, In fact, for a fixed

_{n}*l*this probability is invariant with respect to permutations of the alleles in this sequence [

_{a}*i.e.*, with respect to the way the coalescence events are interspersed (Jenkins and Song 2011)]. There are ways to intersperse the first

*n*−

*l*− 1 coalescence events (Jenkins and Song 2011, Lemma 3.1).

_{a}Now, given , we require the probability that a single mutation event occurs on the correct branch of the tree, an event we denote by Since mutation events occur along the branches as a Poisson process of rate *θ*/2, we have thatas *θ* → 0, where is the total branch length.

Putting all this together,

Summing over *n _{a}*, we findand hence, noting thatsubstituting for the numerator and denominator and letting

*θ*→ 0 yields the given result.□

Henceforward we denote the trinomial coefficient by

*Proof of Lemma 1*. The argument here is very similar to that of the *Proof* of (2). This time a coalescent tree compatible with **v _{n}** and with must exhibit the following sequence of events (Jenkins and Song 2011, Lemma 4.1):

*n*−

_{a}*l*coalescence events of type

_{y}*a*alleles,

*n*−

_{b}*m*coalescence events of type

*b*alleles, and

*n*− 1 coalescence events of type

_{c}*c*alleles, in some interspersed ordering; a mutation event reverting the sole remaining type

*c*allele to a type

*b*;

*l*−

_{y}*l*coalescence events of type

_{o}*a*alleles and

*m*coalescence events of type

*b*alleles; a mutation event reverting the sole remaining type

*b*allele to a type

*a*; and coalescence of the remaining

*l*+ 1 type

_{o}*a*lineages to a most recent common ancestor of the sample.

Here, *l _{y}* is the number of extant type

*a*lineages at the time of the younger mutation event,

*l*is the number of type

_{o}*a*lineages at the time of the older mutation event, and

*m*is the number of type

*b*lineages at the time of the younger mutation event (Figure 3). Arguing as above, any one of the compatible sequences satisfiesand the two mutation events land on the correct pair of branches with probabilityHence, following the reasoning that led to (A1), we findSubstituting for each term in the summand and simplifying leads to (4), and summing over triallelic configurations

**n**yields (6).

The nonnested case is similar. There are possible (Jenkins and Song 2011, Lemma 4.2), each with probabilityThe two mutation events land on the correct pair of branches and are in the correct age order with probability (A2)The additional factor on the right-hand side of (A2) accounts for the fact that, should the two mutations arise during the same epoch (*m* + *l _{y}* + 1 =

*l*+ 1), then only with probability 1/2 is the mutation giving rise to the

_{o}*b*allele the elder one. Finally,which yields (5) after substituting and simplifying. Summing over triallelic configurations

**n**yields (7).□

*Proof of Theorem 2*. Conditional on the event , the single simultaneous mutation event occurs uniformly on the *n* − 1 coalescence nodes of the coalescent tree. Hence, immediately prior to (*i.e.*, more recently than) the coalescence event at which the mutation event occurred, there are *L* lineages ancestral to the present-day sample, with *L* ∼ Uniform{2, … , *n*}. Given *L* = *l*, it is well known (*e.g.*, Griffiths and Tavaré 1998) that the distribution of the number of leaves subtending each of these *l* lineages is uniform on the possible compositions of *n*. Of these, compositions have *n _{b}* leaves and

*n*leaves, respectively, subtending the two lineages coalescing into the node that experiences the simultaneous mutation event. Hence,for 3 ≤

_{c}*l*≤

*n*(and 0 otherwise). Also note thatThis is the probability that the simultaneous mutation event does not occur at the oldest of the coalescence nodes (

*l*= 2), which would lead to a sample containing no copies of the ancestral allele.

Putting all this together we havewhich simplifies to Equation 9.□

## Footnotes

*Communicating editor: J. Wakeley*

- Received October 10, 2013.
- Accepted November 6, 2013.

- Copyright © 2014 by the Genetics Society of America