Notes Lecture on Gene Genealogies
Alan R. Rogers1 All rights reserved. March 22, 2004
1
Department of Anthropology, University of Utah, Salt Lake City, UT 84112
Lecture 6
The Site Frequency Spectrum
6.1 The empirical site frequency spectrum
In a sample of g genes, a polymorphic site may divide the sample into 1 mutant and g 1 nonmutants, into 2 mutants and g 2 non-mutants, and so on. There may be at most g 1 copies of the mutant if the site is to be polymorphic. In many cases we cant tell which allele is the mutant, so category i gets conated with category g i. Such spectra are called folded. I will call a site a singleton if the mutant is present in a single copy, a doubleton if it is present in two copies, and so on.
6.1.1
An unfolded spectrum
Consider the set of DNA sequence data below: 123456 HumanSequence1 AATAGC HumanSequence2 ..AC.. HumanSequence3 .TACT. HumanSequence4 ..ACT. --------------------ChimpSequence1 AAAATC ChimpSequence2 AAAATC There are 4 human sequences and 2 chimpanzee sequences. There are 6 sites of which 4 are polymorphic (segregating) within the human sample. We calculate the empirical spectrum by considering the sites one at a time. Site 1 is xed and therefore does not contribute to the spectrum. Site 2 has both an A and a T within the human sample, but has only an A within the chimpanzee sample. The odds are that the ancestor of humans and chimps had an A at this site, so we can infer that T is the mutant allele. Since there is only one copy of T in the human sample, site 2 is a singleton. So far, our spectrum looks like this: Singletons : 1 25
LECTURE 6. THE SITE FREQUENCY SPECTRUM
26
Doubletons : Tripletons : Site 3 is like site 2. The human sample has a T and 3 As, and the chimp sample has only As. We infer that T is the mutant allele and count this site as another singleton. The spectrum now looks like this: Singletons : 2 Doubletons : Tripletons : Site 4 has an A and 3 Cs, but it appears that A was the ancestral allele. We count this site as a tripleton, so the spectrum becomes Singletons : 2 Doubletons : Tripletons : 1 Site 5 has 2 Gs and 2 Ts. It does not matter which of these is ancestral. Either way, the site is a doubleton. The spectrum becomes Singletons : 2 Doubletons : 1 Tripletons : 1 Site 6 does not contribute to the spectrum. We are done. The empirical spectrum has 2 singletons, 1 doubleton, and 1 tripleton.
6.1.2
A folded spectrum
In the preceding section, the chimpanzee sequences were used at each site to infer which nucleotide was ancestral and which was the mutant. Let us now pretend that we have no chimpanzee sequences and therefore cannot tell the the ancestral allele from the mutant. Instead of counting mutants, we will count the rarest (sometimes called the minor) allele at each site. This time, however, I will omit the invariant sites (1 and 6), which do not contribute to the spectrum. Site 2 The rare allele, T, is present in a single copy, so this site contributes to the singleton category just as it did for the unfolded spectrum. Site 3 Ditto: another singleton Site 4 The rare allele, A, is present in a single copy, so this site is a singleton. Recall that it was a tripleton in the unfolded spectrum. Site 5 Another singleton. The folded spectrum looks like this:
LECTURE 6. THE SITE FREQUENCY SPECTRUM
27
Singletons : 3 Doubletons : 1 The only dierence is that site 4, which was a tripleton in the unfolded spectrum, becomes a singleton in the folded spectrum. In general, the ith category of the folded spectrum contains not only category i of the unfolded spectrum, but also category g i, where g is the number of DNA sequences in the sample.
6.2
The theoretical spectrum
The theoretical site frequency spectrum predicts the values in the empirical spectrum. The only theory that we have is one that assumes a stationary population that mates at random, and a mutational process that conforms to the model of innite sites. I will assume initially that we can tell mutants from ancestral alleles so that our spectrum will be unfolded.
6.2.1
A sites position in the spectrum depends on its position in the gene tree
----A-----| |----B--| | ----------| |----------| | ------C-------------
Consider the following gene tree:
Mutations A and C are singletons, whereas B is a doubleton. A mutation that occurs in the most recent coalescent interval can only be a singleton. One that occurs in the next most recent interval can be either a singleton or a doubleton. One in the interval before that can be a singleton, a doubleton, or a tripleton. And so on.
6.2.2
A tree with two leaves has nothing but singletons
------------------| | |----------| | -------------------
To get a sense of how the process works, it helps to start with a tree with just two leaves:
|-- G generations --|
LECTURE 6. THE SITE FREQUENCY SPECTRUM
28
Since the hazard is 1/G, the mean depth of this tree is G generations and the total length is 2G. We expect 2uG = mutations, all of which will be singletons.
6.2.3
A tree with three leaves has the same number of singletons and half that number of doubletons
Now consider a tree with three leaves: --------| |-------------------| | --------| |-------| | -----------------------------|---G/3--|--------- G --------| Generations If we could look at the spectrum just before the most recent coalescent event, it would look just like that of the tree with two leaves: singletons and no doubletons. At the time of the coalescent event, half of these mutations (the ones on the upper branch) become doubletons. There is no further change in the number of doubletons, so the expected of number doubletons in a 3-leaf gene tree is /2. (We dont need to worry that mutation will turn any of our doubletons back into singletons because, under the innite sites model, mutation never strikes the same site twice.) We start this latest coalescent interval with /2 singletons, but then more singletons are added because of new mutations. How many new mutations should we expect to see? The intervals expected length is G/3 generations (see section 4.3), and it contains 3 lines of descent, so the sum of the branch lengths within this interval is (on average) G generations. We therefore expect uG = /2 new singleton mutations. The number of singletons that is added is precisely equal to the number that was lost. Thus, the new spectrum has singletons and /2 doubletons.
6.2.4
The theoretical spectrum for an arbitrary number of leaves
The argument that I used above gets progressively more tedious as leaves are added. It is better to use a dierent argument. I will skip the details here, but the results look like this [?]: Sample size 2 3 4 5 Theoretical spectrum (singletons, doubletons, . . .) , /2 , /2, /3 , /2, /3, /4 Etcetera
LECTURE 6. THE SITE FREQUENCY SPECTRUM
29
It is remarkable that as we increase sample size, the number of mutants in each category doesnt change. We merely add a new category at the right side of the spectrum. To use the theoretical formula with data, we need to substitute some estimate of . We might use the mean pairwise dierence, , or the estimator S = S
g1 1 i=1 i
where g is the number of DNA sequences in the sample. Either of these estimators might work, since both of them estimate under the stationary neutral model (see the discussion of equation 5.6 on page 22). To choose between these estimators, we need some additional criterion. The sum of the observed spectrum is equal to the number S of segregating sites. It would be useful if the theoretical spectrum summed to the same value. This turns out to be so only if S is used to estimate .
6.2.5
Folded theoretical spectra
When the spectrum is folded, we cannot distinguish category i from category g i. Consequently, the expected number in category i in the folded spectrum is the sum /i and /(g i), the expected numbers in the two corresponding categories in the unfolded spectrum. This works so long as i and g i are not the same number. If they are the same number, then the expected spectrum is simply /i.
6.3
Human site frequency spectra
Figure 6.1 shows all of the human site frequency spectra that I was able to cull from the literature in the year 2000. In each plot, the empirical spectrum is shown as a histogram, and the expected values under neutral evolution with constant population size are shown as bold dots. The top row shows three systems in which there is an excess of singletons, compared with the stationary neutral model. The middle row shows three systems that seem to t the neutral model, and the bottom row shows three systems in which there is a decit of singletons and an excess of sites at intermediate frequency.
6.4
Exercises
EXERCISE 61 Using handouts that I gave out during my rst lecture, (1) use S to estimate , (2) from this value, calculate the number of sites expected in each frequency category, (3) fold the resulting theoretical spectrum by summing values for i and g i. (4) Compare the result with the empirical spectrum in the handout.
LECTURE 6. THE SITE FREQUENCY SPECTRUM
30
MtDNAa,b g = 636 s = 226 Frac. of sites
Yb,c g = 718 s = 20
Xq13.3d g = 69 s = 33
0.0
0.5 0 -globine g = 253 s = 33
1 0.0 Misc Xf g = 10 s = 20
0.5 SNPg g = 147 s = 870
Frac. of sites
0.0
0.5 0.0 LPLh g = 142 s = 88
0.5 0.0 Dystrophini g = 250/860 s = 33 PDHA1j g = 35 s = 24
0.5
Frac. of sites
0.0
p
0.5 0
p
1 0.0
p
0.5
Figure 6.1: Site frequency spectra. The open rectangles in each panel show observed spectra; the bold dots show the spectra expected under the innite sites model with no selection and constant population size. g is the number of chromosomes sampled, s is the number of segregating sites, and p is the frequency of the mutant allele (where ancestral state could be determined) or the rarest allele (where ancestral state is unknown). Sources: a [9], b [6], c [26], d [11], e [5], f [17], g [4], h [3], i [34], j [7]
Appendix A
Answers to Exercises
61 That example was of a sample of g = 10 DNA sequences, which had 15 segregating sites. Thus, we can estimate as
EXERCISE
S =
15
9 1 i=1 i
= 5.3
Ill use the symbols vu and vf to represent the unfolded and folded spectra respectively. The unfolded theoretical spectrum (assuming selective neutrality and constant population size) is vu = [, /2, /3, . . . , /9] The rst entry in this vector is the expected number of singleton sites, the second is the expected number of doubleton sites, and so on. Substituting the estimated value of turns this into vu = [5.30, 2.65, 1.77, 1.33, 1.06, 0.88, 0.76, 0.66, 0.59] The folded spectrum is constructed as follows: vf = [5.30 + 0.59, 2.65 + 0.66, 1.77 + 0.76, 1.33 + 0.88, 1.06] = [5.89, 3.31, 2.53, 2.21, 1.06] Thus, we expect 5.89 sites at which the minor allele is present in 1 copy, 3.31 sites at which it is present in 2 copies, and so on. In the real data (see section 1.3, page 2), we had vf = [6, 2, 2, 5, 0] where the hat indicates that these values refer to data rather than from theory. The theoretical and observed spectra are similar, but certainly not identical. There is no inference that can be drawn from this dierence, because our sample is very small.
77