Genetic Drift Simulation

Stochastic simulation of allele frequencies and heterozygosity under the Wright-Fisher model.

Overview

This repository aims to implement a genetic drift simulation in both Python and Julia. The first one follows the algorithm described by Dr. Laura Figueroa-Corona in the course of Evolution and Population Genetics. On the other hand, Julia's implementation is based on [1], which was provided by Dr. Mashaal Sohail for the Human Genomics course at the B.Sc. in Genomic Sciences.

In both cases, the code simulates genetic drift according to the Wright-Fisher model, in which random sampling leads to allele fixation, allele loss, reduced heterozygosity, and coalescence (see Python's implementation).

Wright-Fisher Model

This models a single evolutive force (genetic drift) with the following assumptions:

• Panmictic population
• Constant population size
• Hermaphroditic individuals
• Discrete generations
• No other evolutive force influences the population

Because the sampling is stochastic, the allele frequency depends only on the population size and behaves as a Markov process due to the lack of memory of past states (the next generation frequencies are influenced solely by the current state).

Repository Structure

.
├── LICENSE
├── Manifest.toml
├── Project.toml
├── README.md
├── results
│   └── images
└── src
    ├── genetic_drift.ipynb
    └── genetic_drift.jl

Running the Simulation

Python's Implementation

This uses a multinomial distribution to randomly sample allele counts because it generalizes the binomial distribution and allows tracking which alleles have been inherited to the next generation. The latter feature was added for a coalescence exercise.In addition, the heterozygosity is computed as follows:

$$ H_e = 2p(1-p) $$

Julia's Implementation

This one uses the traditional binomial distribution as it was not intended for tracking the individual copies each allele inherits to the next generation. In contrast with the Python's version, the heterozygosity $H_e$ is computed as an exponential decay function of time (where time is discrete and corresponds to a given generation) [1]:

$$ H_e(t) = H_0(1-\frac{1}{2N})^{t-1}, \text{ }t > 0 $$

Usage example:

julia --project=. ./src/genetic_drift.jl \
    -i 10 -p 0.3 \
    -g 50 -t 0.42 \
    -o results/images

For more information on the parameters run:

julia --project=. ./src/genetic_drift.jl -h

It is equivalent to:

julia --project=. ./src/genetic_drift.jl --help