Simulates the Wright-Fisher Model of selection and random genetic drift

gistfile1.txt

# data.frame to be filled
wf_df <- data.frame()

# effective population sizes
sizes <- c(5, 10, 100, 500)

# starting allele frequencies
starting_p <- c(.01, .1, .5, .8)

# number of generations
n_gen <- 10

# number of replicates per simulation
n_reps <- 50

# run the simulations
for(N in sizes){
  for(p in starting_p){
    p0 <- p
    for(j in 1:n_gen){
      X <- rbinom(n_reps, 2*N, p)
      p <- X / (2*N)
      rows <- data.frame(replicate = 1:n_reps, 
                         N = rep(N, n_reps), 
                         gen = rep(j, n_reps), 
                         p0 = rep(p0, n_reps), 
                         p = p)
      wf_df <- bind_rows(wf_df, rows)
    }
  }
}

# plot it up!
p <- ggplot(wf_df, 
            aes(x = gen, 
                y = p, 
                group = replicate)) +
  geom_path(alpha = .5) + 
  facet_grid(N ~ p0) + 
  guides(colour=FALSE)
p

gistfile2.txt

#-----------------------------------------------------------
  ## Simulate Genetic Drift (using Wright-Fisher model)
  # http://statisticalrecipes.blogspot.com/2012/02/simulating-genetic-drift.html
  library(ggplot2)
  library(reshape)
  
  # Set up parameters
  N = 2 # number of diploid individuals
  N.chrom = 10*N # number of chromosomes
  p = 0.5
  q = 1-p
  N.gen = 100 # number of generations
  N.sim = 5 # number of simulations
  
  # Simulation
  X = array(0, dim=c(N.gen, N.sim))
  X[1,] = rep(N.chrom * p, N.sim) # initialize number of A1 alleles in first generation
  for(j in 1:N.sim){
    for(i in 2:N.gen){
      X[i,j] = rbinom(1, N.chrom, prob=X[i-1, j]/N.chrom)
    }  
  }
  X = data.frame(X/N.chrom)
  
  # Reshape data and plot the 5 simulations
  sim_data <- melt(X)
  ggplot(sim_data, 
         aes(x = rep(c(1:N.gen), 
                     N.sim), 
             y = value, 
             colour = variable)) +
    geom_line() + 
    ggtitle("Simulations of Genetic Drift") + 
    xlab("Generation") + 
    ylab("Allele Frequency") + 
    ylim(0,1) + 
    labs(colour = "Simulations")

gistfile3.txt

  #---------------------------
  library(learnPopGen)
  genetic.drift(p0=0.5, Ne=20, nrep=10, time=100, show="p", pause=0.1)
  ds <- drift.selection(p0=0.01,Ne=100,w=c(1,0.9,0.8),ngen=200,nrep=5)
  

simulate-Wright-Fisher-Model-of-selection-and-random-genetic-drift.R

# From https://public.wsu.edu/~gomulki/mathgen/materials/wrightfisher.R

# This simualates the Wright-Fisher Model of selection and random genetic drift 
# for an asexual organism with two alleles and no mutation 
# By R. Gomulkiewicz (5 September 2013)
# Updated 8 Sep 2015

s <- 0.01 		# selection coefficient of the mutant type
init.num <- 1	# initial number of the mutant type
N <-10		# population size N
reps <-	15		# number of replicate simulation runs
gens <- 40 		# generations of evolution to simulate for each replicate run
 
#create a matrix with gens+1 rows  and reps columns whose every entry is the initial number of mutants
#below,we will replace entries 2 through gens+1 of each column with simulated numbers of mutants in each replicate run
j=matrix(init.num,gens+1,reps) 

#loops that execute the replicate simulations
#uses simplified form of the post-selection expected frequency pstar = j(1+s)/[j(1+s)+(N-j)] = j(1+s)/(js + N)
for(k in 1:reps){ #this loops over the number of replicate runs
for(i in 2:gens+1)	{  #simulate W-F model starting with init.num mutants
	p.star <- j[i-1,k]*(1+s)/(j[i-1,k]*s+N)  #compute the post-selection expected frequency, given j
	j[i,k]=rbinom(1,N,p.star)  #generate the next j as a single binomial random variable with parameters N and p.star
	}
}

#Show the results as mutant frequencies (i.e., plot each mutant count divided by 2N)
#Display all the replicate runs on a single plot
matplot(0:gens,j/N,type="l",ylab="freq(A)",xlab="generation")