Thursday, August 9, 2012
Measuring population differentiation in R
Biologists often want to be able to measure the degree to which a population is divided into smaller sub-populations. This can be an important thing to quantify, because sub-populations within highly structured populations are, to some extent, genetically distinct from other sub-populations and therefore have their own evolutionary histories (and perhaps futures).
To illustrate this point I've run some simulations. Imagine if we had 5 subpopulations, each with a thousand individuals. In each population we will follow the fate of a locus with two alleles, R and r that have no effect on survival or reproduction and start with frequencies 0.8 and 0.2 respectively (these numbers motivated by this post). In the absence of gene flow between these populations (Panel 1) the frequency of the r allele bounces around due to genetetic drift (evolutionary change, after all, is inevitable). Crucially though, changes in one population can't effect other populations so we end up with substantial among-population differences in allele frequency. In the next two panels, in each generation a proportion of each population's individuals (0.001 and 0.01 respectively) are drawn from the other populations in the simulation. Now that the populations are sharing genes the lines that represent their allele frequencies pull together (that is, the among-population variation is reduced).

One way to quantify the among-population variation displayed in these simulations is to look at the number of heterozygotes you expect to observe across the entire population. The final values for P(r) in the first simulation were {0.33, 0.47. 0.88. 0.10. 0.33} with a mean frequency of 0.42 (so the frequency of the R allele would be 0.58). Knowing our Hardy Weinberg, if we had one big population with two alleles, one being at a frequency of 0.42 we'd expect to get 2pq = 2 * 0.42 * 0.58 = 0.40 heterozygotes. We can call that number HT for expected total heterozygosity. But thats not what we'd actually see in this case. The sub-populations that make up this larger population have their own allele frequencies, when we calculate the expected proportion of heterozygotes for each of these populations by themselves we end up with {0.44, 0.49, 0.21, 0.18, 0.44} for a within-population expected heterozygosity (HS) of 0.35*. This lack of heterozygotes within sub-populations compared with the total population expectation will always arise when genetic drift makes sub-populations distinct from each other. Masatoshi Nei used this pattern to propose a statistic to quantify population divergence called GST, which he defined like this:
Although the problem of GST's dependence on HS has been known for a while, it's taken some time for new statistics that get around this problem to be developed. Philip Hedrick (doi: 10.1554/05-076.1) along with Patrick Meirmans (doi: 10.1111/j.1755-0998.2010.02927.x) introduced G''ST - a version of GST that is corrected for the observed value of HS as well as the number of sub-populations being considered. Meirmans used a similar trick to define φ'ST (doi: 10.1111/j.0014-3820.2006.tb01874.x), another FST analogue that partitions genetic distances into within- and between-population components. Most recently, Lou Joust introduced an entirely separate statistic, D, that directly measures allelic divergence (doi 10.1111/j.1365-294X.2008.03887.x).
The statistical programming language R is becoming increasingly popular among biologists. Although there is a strong suite of tools for performing population genetic analyses in R, code to calculate these "new" measures of population divergence have not been available. My package, mmod, fills this gap. I won't give too many details of the package here, as that's detailed in the paper and the package is will documented. Briefly, mmod has functions to calculate the three statistics described above (and Nei's GST ), as well as pairwise versions of each statistic for every population in a datastet. It also allows users to perform bootstrap and jacknife re-sampling of datasets, the results of which are returned as user-accessable objects which can be examined with any R function (there is also a helper function to easily apply differentiation statistics to bootstrap sample and summarise the results) . The library is on CRAN, so installation is as easy as typing "install.pacakge("mmod")", the source code is up on github. If want to use the package I'd suggest reading the vignette ("mmod-demo") before you dive in.
I'm keen to hear about bugs or feature requests from users, just email them to david.winter@gmail.com
Reference:
Winter, D.J. (in press). MMOD: an R library for the calculation of population differentiation statisticsMolecular Ecology Resources : dx.doi.org/10.1111/j.1755-0998.2012.03174.x
Labels: my research, population genetics, sci-blogs, science
Wednesday, July 11, 2012
You can't ban redheaded sperm
It seems quite a few students have planned their high school and unversity careers in the hope that studying biology meant that could leave maths behind. So, when they are confronted with "p2 + 2pq + q2 = 1" and asked to do something with it, they are unhappy.
That particular formula is for something called the Hardy-Weinberg equilibrium and a significant proportion of students roll their eyes and slump their shoulders when you tell them they are going to need to use it for a problem. They think its arbitrary and irrelevant to anything the least bit important, and what's more it looks a little like it's already solved. So, I'm always looking for ways to convince people that Hardy-Weingberg isn't just simple, but actually intuitive and important. So here's my attempt to explain why knowing a little population genetics is helpful.
You may remember last year a Danish sperm bank had started turning away would-be donors with red hair, since there is little demand for sperm that might contribute to the creation of a readhead. It turns out, if you know a little bit about population genetics you can see that policy will have little effect on the number of red heads the sperm bank helps to bring into the world.
Hair colour is partially controlled by a gene called MC1R. There are different versions of MC1R floating around in human populations, and one of them has a mutation that stops melanin (a dark pigment) passing into hairs as they grow. Geneticists call different versions of a gene "alleles", so we'll call this flavour of MC1R the "red hair allele" and give it the symbol r.As I'm sure you know, you have two copies of most of your genes, one inherited from your mother and the other from your father. Red hair is a recessive trait, which means in order to have red hair both of your copies of MC1R need to be the r allele: if you have one or two copies of the "normal" MC1R allele (which we'll call "R") you have some pigment passing into your hair and it will be another colour. We call the total genetic make-up of a person their "genotype", and their physical characteristics their "phenotype", so here's a table showing the genotypes and the phenotypes we're talking about in this post:
| Genotype |
r/r
|
r/R
|
R/R
|
Phenotype (hair colour)
|
Red Hair
| Not Red Hair | Not Red Hair |
I know there are a lot of technical terms there (Carl Zimmer will not be happy...), but we do need to be precise when we talk about genetics because, strange as it may seem, there isn't a single definition of the word gene. Once you've got your head around the terms, it's all pretty straight forward: you need two copies of the r allele to have red hair. Think what this means for the Danish sperm bank though. Turning away red headed sperm donors doesn't turn away red headed sperm since there will still be "carriers" with only one copy of r (and, thus, non-red hair) donating sperm and half of those sperm will be "red headed sperm".
How big a problem is this likely to be? First we need to work how common the r allele is, and we can use the frequency of redheads to find that. By long tradition, the frequency of a recessive allele is denoted by "q", so, in a population where one quarter of the alleles are r we'd say q = 0.25. We know that in order to have red hair you need both your copies of the MC1R gene to be the r allele and that you inherit each allele separately. When probabilities are independant we can mutiply them, so the chace someone in this population is a redhead is q x q = q2 = 0.06 .Following the same logic, the frequency of the R/R genotype must be the frequency of R squared (by convention, the frequency of a dominant allele is called "p", so that's p2).
Knowing this relationship, we can work backwards and find the frequency of r if we know the proportion of redheads in a population. In most of Northern Europe, about 4% (0.04) of the populaiton are redheads so q2 = 0.04 and q = √0.04 = 0.2. As you can see, red hair genes can be a lot more common than redheads:
It's easy to calculate how the policy would work in populations with more or less redheads:
A lot of rare diseases are caused by recessive alleles. They remain rare for the obvious reason that people with such diseases are unlikely to pass on their genes. But they remain present in populations because, as we've found, once recessive alleles get rare the overwhelming majority of them are found in carriers. In this way, rare recessive alleles are seldom exposed to selection so they stick around for a long time.*
Because disease causing genes stick around in populations, there is a pretty good chance that you carry a few alleles that would cause a debilitating disease in someone who had two copies of them. The same applies to anyone you might be hoping to have children with. Thankfully, its very unlikely that your prospective mate with have disease-causing alleles for the same genes that you do. That is, as long as you look beyond the family tree when you look for a mate. If you have a child with someone who is closely related to you, you will have each inherited some of your genes from the same source, which increases the chances you share disease alleles.
*In fact, they often reach a point called "mutation-selection balance" in which the frequency of the allele remains static because new mutations re-create the allele as quickly as selection removes it. JBS Haldane was the first person to notice this, and he used his theory to create a very accurate estimate of the human mutation rate well before we knew what genes were made of!
Labels: genetics, population genetics, sci-blogs, science, teaching
Tuesday, October 27, 2009
Some updates
- My support for the Bellbird in the annual Bird of Year poll was not enough for it raise above a mid-table finish. In a fit of creativity kiwis voted the Kiwi as their favourite bird with rifleman and kea filling the minor placings.
- That nutty paper about insect metamorphosis was held up at PNAS while the editorial board looked into Lyn Margulis' own claim that she went through seven reviewers to find two that would sign off on it. Even so, the paper will make it to a dead-wood issue of PNAS along with a letter from Gonzalo Giribet who told a Nature reporter the idea was "the most stupid thing that has ever been proposed"- should be fun.
- Finally, if you thought that the idea of their being some genetic differences between human races wasn't a controversial idea then you might check out two threads (1, 2) that started when Larry Moran talked about the same article at The Sandwalk.
Labels: bellbird, evolution, human genome, Korimako, population genetics, sci-blogs, science
Friday, July 31, 2009
My mutant mitochondria and life on earth
You may know Damian from conversations at Ken's and various places that approach similar questions. You should probably go read his blog because it's interesting and, frankly, he's desperate for the traffic. A while back Damian was thinking about mitochondria and what they can tell us about people. Like any good scientific pedant I turned up and picked holes in specifics of what Damian had to say but provided nothing comparable to the message he was trying to provide. So now, only three months after Damian provided his thoughts on mitchondria, here's a few of thoughts on my own mitochondria, the dark secret each one holds and why that secret is the key to understanding just about everything in biology.
But, before I out myself to the world I should perhaps talk about just what a mitchondrion is. At it's most basic level all life is a very special set of chemistry, our cells are chemical factories making the materials from which we are built and running the reactions that keep us ticking over. 1 The cells of complex organisms like you, me and slime molds can run thousands of different chemical reactions at the same time without the products of each reaction interfering with each other because we have structures called organelles (yes, small organs) that effectively wall off one bag of chemicals from the rest of the cell. Mitochondria are the organelles that make energy for the rest of the cell's activities. Which really ought to be reason enough to celebrate them in a blog post, but as an evolutionary biologist I have another reason to venerate mitochondria. In animal and fungal cells mitochondria are the only structure outside of the the nucleus that contain their own DNA and in most animals that DNA (which I'm going call mtDNA from now on) is passed from mothers, and not fathers, to offspring (analogous to the way surnames pass from fathers and not mothers). To me a mitochondrion is not only a source of cellular energy but an independent witness to the evolutionary history of its host.
So, what can this extra set of genes tell us? When I was an undergrad I spat in a falcon tube, extracted DNA from the check cells that I had floating around in my mouth and amplified some of my own mtDNA then sent the sample away to have its sequence of chemical bases read and converted into a sequence I can read:
AAGGAACAATCAGAGAAAAAGTCTTTAACTCCACCATTAGCACCCAAAGCTAAGATTCTAATTTAAACTATTCTCTGTTCTTTCATGGGGAAGCAGATTTGGGTACCACCC
AAGTATTGACTCAGCCCATCAACAACCGCTATGTATTTCGTACATTACTGCCAGCCACCATGAATATTGTACGGTACAATAAATACTTGACCACCTGTATACATAAAAACC
CAATCCACATCAAAACCCCCTCCACATGCTTACAAGCAAGTACAGCAATCAACCCTCAACTATCACACATCAACTGCAACTCCAAAGCCACCCCTCACCCATTAGGATACC
AAACAAACCTACCCACCCTTAACAGTACATAGCACATAAAGCCATTTACCGTACATAGCAACATTACAGTCAAATCCCTTCTCGTCCCCATGGATGACCCCCCTCAGATAG
GGGTCCCTTGACCACCATCCTCCGTGAAATCAATRATCCCGCCACAAGAGTGCTACTCTCCTCGCTCCGGGCCCATAACACTTGGGGGTACGCTAAAAGTGAACTGTATCC
GACAATCTGGTTCCTACTTCAGGGGCCATAAAGCCTAAATAGACCCACAACGTTCCCCCTTAAATAAGACCATCACGATGGATCACAGGTCTATCACCCTATTAAACCACT
CACGGGGAGCTCTCCATGCATTTGGGTATTTCGTACCTGGAGGGGGTATGCACGCGGATAGCATTGCGAGACGCTGGAGCCGGAG
There it is then, highlighted in blood red, my dark secret. The prosaic explanation for that bright red 'R' amongst the 'A's 'T's' C's and 'G's is that the sequencer couldn't decide if I had a 'A' or a 'G' in that position. On closer inspection it turned out I had both - some of my mitochondria have 'A's, some 'G's. Somewhere in the line of mothers and daughters that lead to the source of my mitochondria, or in my Mum, or possibly even sometime during my early development an error was introduced to one mitochondrion's DNA and over time that error was copied so many times that by now it's present in about half of my mitochondria. That is to say, I am a mutant.
Perhaps I'm being overly dramatic. The region of mtDNA that I sequenced doesn't make a protein, it's not quite junk DNA but the particular DNA base at the highlighted red position probably has no discernible effect on the way my cells work. Our typical conception of mutation is drawn from the tragic effects of those relatively rare mutations, induced in our bodies or passed on through germ cells, that lead to diseases (or, in movies to super powers). In fact, we are, each of us mutants. DNA replication is not perfect, we are born with about 6 or 7 new mutations [err, I was out by a factor of 20 there, more like 100-200 mutations...]in our nuclear genome and the mitochondrial genome mutates much more quickly than that. One of the revelations of evolution biology in the molecular age has been the realisation that most mutations are like the one I've revealed to world - of little or no effect. They occur in regions of the genome that don't have genes or if they are in genes they don't alter that gene's product or even they do alter the product the difference is so slight the end result is undetectable. Only a few mutations have devastating effects, an even smaller minority actually make live easier for their host.
At first glance the genuinely silent majority represented by 'neutral mutations' might seem a more boring topic than their deleterious and advantageous counterparts ( the molecular basis of, respectively, disease and natural selection). In fact, because neutral mutations behave in predictable ways within populations the development of a neutral theory of molecular evolution has allowed us to test a lot of ideas in evolutionary biology. Let's look at an example. In each generation every individual has a small chance of having a mutation in at each of their DNA bases, we call this the mutation rate (µ). This means in a population new mutations arise at a given base at the individual mutation rate multiplied by the population size (N). What happens to these mutations once they get into populations? If they are selectively neutral they won't effect their host's chance of reproducing so it will come down to chance - if the mutation arises in a particularly fecund individual there will be lots of copies of the mutation in the next generation, if our mutant is struck by lightning before they reproduce the mutation will go extinct. As long as the mutation survives its frequency in subsequent generations will bounce up and down with no particular direction but in the long run finite population sizes mean there are only two fates for mutations, they go extinct or they completely take over the population (become fixed in population genetics parlance). We can calculate the probability that a new mutation becomes fixed rather than lost, it's 1/N. Where did I pull that from? A population genetics model called the coalescent might help to explain why this is true.
Take a very small population, six individuals:
You might be wondering what kind of trick I'm up to here, offspring don't choose their parents and even if we are doing it at random aren't we going backwards? Yes, but that's fine. We are not describing populations, we're making a model to help us understand how they work. Saying that each child has a 1/6th chance of we're being the offspring of each parent (and deciding that with a die) is exactly the same thing as saying each potential parent has a 1/6th chance of being the parent of each member of the next generation. So, even if offspring choosing their parents is biological nonsense it's a useful way of understanding real biology. With that out of the way lets look at our parental generation, the 4th and 5th individuals didn't reproduce. I'm sure the led rich and full lives and contributed to society in many interesting ways but at the moment we're interested in how we ended up with our current population, so lets discard them then simulate a few more generations, ignoring individuals that didn't contribute to the last generation:
Ha! Our lineages have coalesced. In fact, I shouldn't be surprised because every lineage for every gene in a population will eventually coalesce at a single common ancestor (called the most recent common ancestor or MRCA). There is a MRCA for all human mitochondrial lineages, kiwi evolutionary biologist Allan Wilson named her mitochondiral eve which I'm sure he thought was very clever at the time but has since led many people to that mtEve was, in every way, like her biblical counterpart. You can see from our little simulation that she need not have been the only human on earth and, in fact, there were other people alive at the same time as mtEv who also have modern descendants (because eve is only the MRCA of the matrilineal line). Moreover, in large populations the title mtEve can only be bestowed on someone who has been dead for many, many generations and in subsequent generations as lineages die out, the title will pass on to someone else.
Let's leave Eve for the time being and return our focus to our population of circles. Think about what the inevitability of coalescence means for new mutations. Each individual in generation 5 will either be the ancestor of all the individuals in the present generation or none of them. Under the neutral model those eventualities are decided at random so the chance that any gene (including a new mutant) becomes the MRCA as some later stage really is one over the population size - 1/N. Now we have a term that describes the rate at which mutations arise and another that tells us how likely they are to be fixed once they have. From here it one tiny step to work out the rate at which new mutations will be fixed (k) in a real population:
k = (µ x N) x 1/NAnd we can get rid of the population size in a puff on third form algebra:
k = (µ x N) / N = µThat is, new neutral mutations are fixed in populations at a constant rate equal to the population mutation rate (and not effected by population size). Of course the individual mutation is probabilistic, so the 'constant' rate is in fact the cumulative adition of discrete events meaning the actual number of mutations fixed in a given population at a given time will vary around the expected value. Having an idea of how a real population behaves in the absence of selection provides evolutionary geneticists with a null hypothesis to test for selection. Let's look at some real data from a the human and chimp genome projects. You may already know that most of the DNA in mammalian genomes is not actually made into a gene product - even in the regions of genome that actually make proteins are flanked by "untranslated regions" (UTRs) that are required for the gene to work but don't code for part of a protein2. Ryuichi Sakate and colleagues compared the rate at which mutations have been fixed in the coding regions and the UTRs of genes in either the human or chimp genome since we parted ways with our cousins
As you can probably see mutations are fixed signifcantally more frequently in the UTRs of the genes than they are in the coding sequence (an effect that is slightly masked by the fact that translated regions are smaller than coding sequences so are more likely to have no mutations at all). Remember, that our neutral model tells us that mutations are fixed at a rate equal to the individual mutation rate. There is no reason to believe that the mutation rate of UTRs is any greater than that of the coding region that they are associated with so the the fact less of less those mutations are getting fixed must be down be differences in the survival of those mutants - natural selection. In fact, in this case we have good evidence for 'purifying selection' weeding out deleterious mutations in the coding regions of genes with a more permissive selection regime in the UTRs since mutations ocuring here can't effect the sequence of the protein that the gene produces.
Even in the coding region only a subset of mutations can change the protein sequence so we can further classify the mutations in Sakate et. al's study into 'silent mutations' that don't change the protein sequence (dS) and 'non-silent mutations' (dN) which do change the protein. Due to the nature of the genetic code we'd expect silent and non-silent mutations occur at about equal rates. When you look at the inset of the graph above you find that, for these genes at least, silent mutations are much more likely to be fixed than non silent ones. Again, this is good evidence that we have purifying selection, most of the non-silent mutations that occur in the coding regionsare being selected against. This probably shouldn't come as a great surprise, after all our genes are the results of round upon round of natural selection - they are already pretty good at what they do so any change is likely to make the protein worse 3. What's important is that the very simple things we deduced about mutations above allows us to understand how natural selection has worked in the human genome.Or am I pulling a fast one? The analysis above depends on the idea that humans and chimps once shared a common ancestor and as well all know that is a wildly speculative hypothesis based on tenuous extrapolation from a few partial fossils and a lot of wishfull thinking from people whose faith in atheism is so great they need to banish a creative spirit from their lives. Or perhaps not. Think back to our population of circles, the key message from that simulation was that within a population mutations are fixed at a constant rate. The generation of new species, which we call speciation, is the splitting of a single population into two distinct lineages that no longer share genes.
During and after speciation, each new species will start accruing new mutations independently and, for neutral mutations, at a constant rate. As the new species continue on their new trajectories their DNA sequences will become progressively more distinct from each other as they independently fix mutations. We can use this information to discover relationships between species. Let's go the The Big Gene Database and get some sequences that match a subsectin of my mitochondrial DNA:
me CATTACAGTCAAATCCCTTCTCGTCCCCATGGATGACCCCCCTCAGATAG chimp CATTACAGTCAAATCCATCCTCGCCCCCACGGATGACCCCCCTCAGATAG gibbon CATCCCAGTTAAATC-ATCCTCGTCCCCACGGATGCCCCCCCTCAGATGG monkey CATATTCATTAAATA-ATCCTCTTCACCACGGATGCCCCCCCTCACTTAG
Now, calculate the proprtion of sites which are different between each sequence:
| me | ||||
| 0.082 | chimp | |||
| 0.163 | 0.122 | gibbon | ||
| 0.306 | 0.265 | 0.224 | monkey |
So, the smallest difference is between me and the chimp. If you keep doing this process, but now recording thepercentage difference between either my sequence or the Chimp one against the Gibbon and Old World Monkey sequence you find that the Gibbon is more similar to the human-chimp group than it is to the Old World Monkey sequence. We can represent these relationships with a a tree (which might make this blog's logo make sense to you):
This result is the antidote to people that make the claim evolutionary biologists are simple comparing 'similarity' when we estimate the relationships between species. Here we have used a simple model of the way in which sequences change and applying that to the sequences we have to work out the relationship which most easily predicts the differences between those sequences. In actuality our model is overly simplistic (for instance certain DNA changes are more likely to happen than others so we can't use raw percentage difference between sequences) and using only 50bp of DNA to estimate relationships is unlikely to get you published in Nature.The best evidence that the phylogeny produced above is the eight one is the fact independently lines of evidence (unlinked genes, morphology, geography, behaviour) support it.
There is one last thing that we can work out from our DNA sequences and what we know about population genetics. If we can put a date on one of the branching points on the tree we can actually estimate the time for all the other splits. As is happens we know that the oldest Old World Monkey fossil to be discovered is around 15 million years old, meaning the latest that first branch that splits the apes (me, the Chimp and the Gibbon) from the monkey sequence could have happened is 15 millions years ago. The average difference between the monkey sequence an a ape one is 0.265 - which divided by 15 million years gives us substitution rate of around 0.018 bases per million years. If we apply this rate to observed distance between human and chimp sequences (0.082 changes per base /0.018 base changes per million years) you get an estimate of the age of the split at around 4.5 million years. When you do the same analysis on more than 50 bases of DNA you normally end up with a bunch something more like 5-7 million years.
Sadly my mutations is, from an evolutionary perspective, effectively dead. Males don't pass on their mitochondria so that mutation will die with me. Still, as I said each carry our own set of private mutations - wouldn't it be nice to think one day a variant that started as something unique to you made it so far through humanity that it could be used to keep track of the way natural selection continues to act on our species and even to help us find our place in the biosphere?
1. Creationists, this a literary device known as "metaphor", please move along2. There are also regions called introns that are spliced out of genes before they moce off to be translated, but we'll focus on the UTRs here
3. There are genes, even in this dataset, for which there is evidence for positive selection - mutations being fixed more quickly than you'd expect under neutrality - which is the basis of adaptatoin to local environmental conditions. But in the main most selection is purifying.
Labels: evolution, might interest someone, mitochondria, mutations, population genetics, sci-blogs

