Coalescent Theory (1) -- genetic events

02 Apr 2015

"All models are wrong. We make tentative assumptions about the real world which we know are false but which we believe may be useful." -- Box

A genealogy is completely summarized by the entire topology (who relates to whom) and the length of each branch (intervals between two subsequent convergence/divergence events). The coalescent approach generates the genealogy backwards, instead of forwards, for a sample of sequences (rather than the entire population). It traces the ancestral lineages, which are the series of genetic ancestors of the samples at a locus, back through time. The history of a sample of size n comprises n − 1 coalescent events. Each coalescent event decreases the number of ancestral lineages by one. This takes the sample from the present day when there are n lineages through a series of steps in which the number of lineages decreases from n to n − 1, then from n − 1 to n − 2, etc., then finally from two to one. At each coalescent event, two of the lineages fuse into one common-ancestral lineage. The result is a bifurcating tree. This approach introduces computational and analytical convenience because the history of the entire population includes sequences that are extinct or we have not sampled.

Gustave Malécot (in the 1940’s) introduced the idea of following a pair of gene copies back to their common ancestor and the notion of identity by descent (IDB): If we pick two genes from a population whose evolution is described by a Wright-Fisher model (discussed below), how long ago on average did the two genes share their most recent common ancestor (MRCA)? Genealogical approaches to samples larger than two appeared in response to the first direct measurements of molecular variation (Harris 1966; Lewontin and Hubby 1966). Since then there have been some seminal work:

In the rest of this post, I will briefly introduce the coalescent approach, including its primary components and extensions. The ultimate goal is to see how this approach can help us understand patterns of variation in DNA sequences. To this end, we will first look into the genealogical process describing how a sample of DNA sequences are related through their common ancestry (i.e., how they coalescent to a single DNA sequence), assuming the process follow the Wright-Fisher model. Then we will relax the model’s constrained assumptions by introducing other important processes. For example, the second process to be incorporated is how different alleles as mutations arise over time. These mutations are the variations that we are interested in. Recombination (cross-over or gene conversion) is another key process that breaks down associations among single sites of a chromosome and allow us to locate interesting variations in preferably shorter regions (see linkage disequalibrium).

In addition to three primary processes, several other genetic or population events will be introduced as important extensions.</b> Together these events (basics and extensions) form a stochastic coalescent process that produces not only a set of simulated DNA sequences but also their genealogical relationship and evolutionary histories. This process is assumed to have so-called Markov property -- the state of a specific next step in the process only depends on the present state (i.e., memory-less). In genetics it is natural to assumes that the probability of a future event depends only on the current situation. When time is continuously measured (as a non-negative real value), the above process is called a continuous-time Markov chain/process and is closely related with the exponential distribution (the only continuous distribution that is memory-less). If a random variable U follows an exponential distribution, denoted U ~ Exp(λ), then Prob(U <= t) = 1 - exp(-λ*t), where λ > 0 is the parameter of the distribution called the rate parameter. Let X1, ..., Xn be independent exponentially distributed random variables with rate parameters λ1, ..., λn.

A property of exponential distributions is that min{X1, …, Xn} is also exponentially distributed, with parameter λ = λ1 + … + λn. The index of the variable which achieves the minimum is distributed according to the law Prob(Xk=min{X1, …, Xn}) = λk/(λ1 + … + λn). As will be more clear below, this property is important for constructing a unified model incorporating various evolutionary events that are assumed to have independent exponentially distributed waiting times. These events compete to occur during a single coalescent process. The time interval between two subsequent same-type events has an exponential distribution parameterized by the rate of this type of events (scaled on population size N). The time interval between two subsequent events (not necessarily the same type) is the minimum of all exponential distributions, which itself is also an exponential distribution with a parameter given by summation. For more information, please refer to the resources provided at the end of this post.

The Wright-Fisher model

The Wright-Fisher model provides a dynamic description of the evolution of an idealized population and the transmission of genes from one generation to the next. This basic model of reproduction makes the following assumptions:

A reproduction process following the Wright-Fisher model results in a decay of genetic variation. Since the population is finite in size and reproduction is a random process, some individuals may not contribute any offspring to the next generation. This random loss of genetic lineages forward in time is called genetic drift, which reduces the diversity of the population diversity. One measure of population diversity is heterozygosity, defined as the probability that two genes chosen at random from the population have different alleles. Alleles are different versions of the genetic information encoded at a location in the genome of an organism (aka, genetic locus). A common example of genetic locus is the sequence of nucleotides that makes up a gene. Thus, two sequences of the same gene are different alleles if they are not identical. Assuming a gene has two allelic states (denoted A and a), genetic drift eventually leads to either A or a being lost from the population. When this happens, the surviving allele is said to be fixed in the population. The effect of genetic drift is compensated by mutation, a process by which the allelic state of a gene occasionally changes from one to another (e.g., from A to a).

Coalescent Events

Parameters (as denoted in the literature/as the input of Cosi or CoJava if applicable)

Algorithm 1

  1. Start with k = n sequences
  2. Simulate the waiting time to next coalescent event, T(k) ~ exp((k(k-1))/2), that is, Prob(T(k) <= t) = 1 - exp(-λ*t) where λ = k(k-1)/2
  3. Choose a random pair (i, j) with 1 ≤ i < j ≤ k uniformly among k(k-1)/2 possible pairs
  4. Merge i and j into one gene and decrease the sample size by one, k--
  5. If k > 1 go to Step 2, otherwise stop

Note: Step 2 utilizes the property that when n is much smaller than N, the probability of a coalescence event in a given generation with k sequences (i.e., for k genes to have k-1 ancestors in the previous generation) is approximately k(k-1)/(4N). Thus, the amount of waiting time (measured in 2N-generation units) during which there are k lineages, T(k), has approximately an exponential distribution with mean k(k-1)/2. The mean of an exponential distribution is equal to the decay rate of the same distribution (i.e., λ). We generate a random number with an exponential distribution (say x) by first generating a uniform random number, u, in [0,1), and then calculating x by: x = ln(1-u)/(−λ), where λ is the rate parameter of the exponential distribution.

Primary CoJava functions
Update the ARG: /populationStructures/demography.java/coalesceByName()

Mutation events

mutation

We consider strictly neutral mutations that will not affect an individual's fitness (the individual's ability to survive and to produce offspring). Such mutations should not affect the simulated genealogies, because they have no effect on the number of offspring or individuals' tendency to migrate. This has two consequences.The first consequence is an efficient computer algorithm, in which the coalescent process is modeled by separating the neutral mutation process from the the genealogical process. We can first generate the random genealogy of the individuals backward in time, and then superimpose mutations forward in time. The second consequence is that we can choose from various mutation models (e.g., infinite-allele, infinite-site, or finite-site model) without influencing the statistical properties of resultant genealogies.

The infinite sites model assumes that a mutation will always happen in a new site/position (thus "infinite sites"), so all mutations are distinguishable (i.e., no recurrent or reverse mutations). Moreover, there will always be one or two states in a position, never more, because each position mutates at most once. Thus, alleles are often labelled as 0 or 1(and a sequence as a string of zeros and ones) regardless their specific meanings. The infinite sites model can be interpreted as describing the evolution of very long DNA sequences with low mutation rate at each position. In contrast, the finite sites model admits that a DNA sequence has a fixed length. Both models assume that positions evolve (via mutation) independently of each other, i.e., a mutation in one position does not influence the chance of a mutation in another position. In principle a mutation model should also describe chromosome-level mutations (e.g., insertions, deletions, etc.), but these events occur so rarely that they normally can be ignored. For the infinite sites model, we can set up an overall rate to comprise both the rates of nucleotide-level and chromosome-level mutations. For the finite sites model, however, usually only nucleotide substitutions are modeled.

Parameters

Algorithm 2

  1. Execute Algorithm 1 to simulate the genealogy of n sequences
  2. For each branch draw a number, Mt, from a Poisson distribution with intensity tθ/2, where t is the length of the branch. Later there will be Mt mutation events added to this branch.
  3. Starting at the root, move forward in time and modify the sequences produced in Step 1. For infinite sites models, add Mt mutations to the descendant sequence of each branch; the position of a mutation is randomly chosen along the sequence.

Primary CoJava classes or functions /geneticEvent/mutations.java /coalSimulator/sim.java/simMutate() By default, Cosi and CoJava is a finite-sites simulation in that mutations occur at discrete sites and if multiple mutations occur at a single site, only the first one is retained. However, setting the parameter “infinite_sites” to yes (see the exemplary param file in this post) converts the output positions to floating point, with all mutations retained. Moreover, both programs allow users to set up the (fixed) number of mutations.

Recombination events

When diploid individuals reproduce, there are two parents, each of which contributes one of its two homologous chromosomes, or a combination of both when they undergo recombination. As shown in the figure below, an individual obtains a chromosome from a parent and the chromosome has two different ancestor sequences, each from one grandparent.
cross-over
As in the standard WF model, we ignore the existence of individuals and focus on chromosomes. Since a recombination event splits the genetic material of a sampled sequence onto two different ancestors (but each single point on the sequence has exactly one thread connecting all its ancestors), it is formulated as opposite and competing to an coalescent event, which combines two sample sequences into one ancestor. The following figure shows the two types of events in a coalescent process.

Parameters

  • r: recombination rate as the probability of a recombination event between two adjacent nucleotide sites per generation (between a parent and a child). Analogous to the definition of the scaled mutation rate, we often define the scaled recombination rate as ρ = 4Nr (population size N, not the sample size n).

    • Algorithm 3: a WF model with recombination

    1. Start with k = n genes.
    2. For k sequences with ancestral material, draw a random number from the exponential distribution with parameter k(k-1)/2 + kρ/2. This is the time to the next event.
    3. With probability (k − 1)/(k − 1 + ρ) the event is a coalescence event, otherwise it is a recombination event.
    4. If it is recombination, draw a random sequence and a random point on the sequence. Create an ancestor sequence with the ancestral material to the left of the chosen point and a second ancestor with the ancestral material to the right of the recombination point. Increase the number of ancestral sequences k by one and go to 1.
    5. If it is a coalescence event choose two sequences among ancestral sequences at random and merge them into one sequence inheriting the ancestral material to both of the sequences. Decrease k by one. If k = 1 end the process, otherwise go to 1.

    Note Assuming lineages recombine independently, the total rate of recombination when there are k lineages is kρ/2. Thus, given only coalescence and recombination events, lineages increase at a linear rate and decreases at a quadratic rate (the total rate of coalescence whenever there are k lineages is k(k − 1)/2). As a result, the number of lineages is guaranteed to stay finite and will eventually hit one. Notably, we also assume that a sequence is unlikely to get involved in a recombination event and a coalescent event at the same time (as illustrated below), because their co-occurrence has a relatively low probability of 1/2N * r = ρ/(8N^2) when r = ρ/(4N), which is negligible for large N and fixed ρ. negligible events
    Actually we often assume away the co-occurrence of different types of events (recall that they are considered as competing events), given that population size (N) is big enough. (Note: there can be zero event at a specific generation). Primary CoJava functions
    Find the recombination location on the sequence: geneticEvents/RecombWorker.java/recombExecute() Update the ARG: /populationStructures/demography.java/recombineByIndex() Create ancestral sequences: /populationStructures/nodeWorker.java/nodeRecombine() Cosi and CoJava allow user to define the recombination rate, which can be to vary according to a predefined recombination file (see this post), which defines a genetic map or hotspots along the genome

    Gene conversion events

    Theoretically, gene conversion is an event in which a portion of the sequence of one chromosome is altered to form a copy of another homologous sequence. It occurs as a result of mismatch repair of double-strand breaks during recombination (and thus is common at recombination hotspots), in which genetic information is copied from an intact sequence to the region of recombination that contains a double-strand break. Despite its complicated and not fully understood biological mechanism, gene conversion is relatively easy to model regarding its consequence: Given two DNA sequences (they are homologous from a genetic point of view), the ‘acceptor’ sequence is replaced, wholly or partly, by a sequence copied from the ‘donor’, whereas the sequence of the donor remains unaltered. Thus, we can model gene conversion as two very close crossover points (using Algorithm 3) even though this rarely is the mechanical way it occurs.

    The figure below sketches the difference between cross-over and gene conversion in backward simulation. Note that there may only be one break point within the sequence if the tract extends beyond the end of the sequence, or if a tract initiates outside the sequence but ends within. In these cases, gene conversion will be indistinguishable from cross-over.
    twoRecombType

    Parameters

  • | geneConversionRate: an input parameter in CoJava
  • 1/Q | gcMeanTractLength: In practice, the exact length of gene-conversion tract (i.e., the portion of the ‘acceptor’ sequence copied from the ‘donor’) usually cannot be precisely known. Empirical evidence show that the tract length (in nucleotides or bp) is drawn from a geometric distribution. In the infinite sites approximation the tract length becomes exponentially distributed with intensity Q such that 1/Q is the mean length of the gene conversion tract.
  • | gcRate: calculated by geneConversionRate * (length of sequence + length of gene conversion tract)

    • Algorithm 4 When a gene conversion event occurs, the first break point is chosen uniformly on the sequence, and the second break point is chosen a distance away from the first as determined by a random number from the tract length distribution.

    Primary CoJava functions
    Find conversion locations on the sequence: geneticEvents/gc.java/execute() Update the ARG: /populationStructures/demography.java/gcByIndex() Create ancestral sequences: /populationStructures/nodeWorker.java/nodeRecombine()

    Migration: User-defined matrix Mating system: Random Mating Fecundity: Random Distribution Life cycle: Not applicable Population growth: Exponential growth Events allowed: change in population size, colonization or extinction (change in number of populations), migration matrix Selection: Single biallelic site

    cannot distinguish…since they have the same total scaled mutation rate. For example, (1) simulate 1,000 sequences from a population of 20,000. Each sequence consists of 2,000 independent loci each of lkb long. Mutation rate is 10-8 per base pair per generation, (2) simulate 1,000 sequences from a population of 1,000. Each sequence consists of 40,000 independent loci each of lkb long with the same mutation rate in (1).

    References