inst/doc/popgen_drift.R

In evolved: Open Software for Teaching Evolutionary Biology at Multiple Scales Through Virtual Inquiries

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
options(rmarkdown.html_vignette.check_title = FALSE)

## -----------------------------------------------------------------------------
library(evolved)
library(plot3D)

# Let's also store our par() configs so 
# we can restore them whenever we change it in this tutorial
oldpar <- par(no.readonly = TRUE)  

## -----------------------------------------------------------------------------
N <- 32 #population size
n_alleles <- 2*N 
p_gen0 <- 0.25 #Frequency of allele A1 in the first gen
p_gen1 <- rbinom(1, n_alleles, p_gen0) / n_alleles
p_gen1

## -----------------------------------------------------------------------------
p_gen2 <- rbinom(1, n_alleles, p_gen1) / n_alleles
p_gen2

## -----------------------------------------------------------------------------
p_gen3 <- rbinom(1, n_alleles, p_gen2) / n_alleles
p_gen3

## -----------------------------------------------------------------------------
p_gen4 <- rbinom(1, n_alleles, p_gen3) / n_alleles
p_gen4

## -----------------------------------------------------------------------------
p_gen5 <- rbinom(1, n_alleles, p_gen4) / n_alleles
p_gen5

## -----------------------------------------------------------------------------
generations <- seq(from = 0, to = 5, by = 1)
p_through_time <- c(p_gen0, p_gen1, p_gen2, p_gen3, p_gen4, p_gen5)
plot(generations, p_through_time, type="l", lwd = 2, col = "darkorchid3",
     ylab = "p", xlab = "generations", las = 1)

## -----------------------------------------------------------------------------
#We will simulate five flasks through 10 generations:
WFDriftSim(Ne = 32, n.gen = 10, p0 = 0.5, n.sim = 5)

## -----------------------------------------------------------------------------
# R will not return anything in the console when you run this, but 
# after running you should have this function in your R session
# as an object. Note that if you name another object as "isFixed"
# the function will be overwritten.

isFixed <- function(p, tol = 0.00000000001){
	
	if (p <= tol | p >= (1 - tol)){
		return(TRUE)
	} else{
		return(FALSE)
	}

}

# Assuming your current allele frequency is stored in the p_t object, the
# line below uses the function created above to test if p_t is equal to 
# zero OR equal to one

## -----------------------------------------------------------------------------
# The function will return `TRUE` if p_t is equal to zero or one.

## -----------------------------------------------------------------------------
WFDriftSim(Ne = 32, n.gen = 10, p0 = 0.5, n.sim = 100)

## ---- echo=F------------------------------------------------------------------
knitr::include_graphics('Buri_1956_fig6_wout_axes.png')

## -----------------------------------------------------------------------------
data = WFDriftSim(Ne = 8, n.gen = 50, p0 = 0.5, n.sim = 50,
                  print.data = T, plot.type = "none")

#then, if we want to see the 2nd to 5th generations of 
# the simulations number 14, 17 and 18, we type:
sims <- c(14,17,18)
gens <- 3:6
subset_data <- data[sims, gens]
subset_data

#and we can also plot those results:

#first opening an empty plot:
plot(NA, ylab = "Alleleic frequency", xlab = "Generation",
     xlim = c(2, 5), ylim = c(0, 1))

#creating a set of colors to paint our lines
cols <- rainbow(n = nrow(subset_data))

#then we finally add the lines to our plot:
for (i in 1:nrow(subset_data)) {
  lines(x = gens - 1, y = subset_data[i,], col = cols[i])
}

## ---- echo=FALSE--------------------------------------------------------------
pops <- c(rep(100000, times= 5), 50, rep(1000, times= 4))
CS <- log(pops, base=10)
time <- 1:10
plot(x=time, y=CS, xaxt="n", yaxt="n", ylim=c(1,5), ylab="Census size (individuals)", xlab="Time since study started (generations)")
lines(x=time, y=CS)

#x ticks:
xtick<-seq(1, 10, by=1)
axis(side=1, at=xtick, labels = TRUE)

#y ticks:
ytick<-c(1,log(50, 10),3,5)
axis(side=2, at=ytick, labels = c("1", "50", "1000", "100000"))

## ---- fig.width=10------------------------------------------------------------
# A simple genetic drift simulator as a Markov process
N       <- 16
popsize <- 2*N

# Now we define a transition matrix, such that:
# Each element of the transition matrix is a transition probability
#   Specifically: element Pi,j is the probability of going from
#   j alleles to i alleles in a single sampling step 
#  

tmat <- matrix(0, nrow = popsize + 1, ncol = popsize + 1) # along rows: current allele count; along cols: future allele count

# Fill in elements with probabilities from the binomial distribution
# If you have j alleles in generation t-1, then this defines the probability
# of sampling (0,1....2N) alleles in the next generation, under a binomial 
# sampling process:

# Note also that first column and row are absorbing states of 0

# Here is the vector of frequencies:
probvec <- (1:popsize) / popsize  

# What is the probability of getting to 0 allele copies, given 
# that you have 1? This is P01, and the transition probability 
# would be computed using 1/2N (=1/32) as the sampling probability 
# for the allele. 
# Here, this is computed using the function dbinom, which is a shortcut 
#   to the analytical binomial probability:
dbinom(1, popsize, prob = 1/32)

# Now we will set up the full matrix by iterating this over rows:

for (ii in 1:nrow(tmat)){
	tmat[ii, 2:ncol(tmat)] <- dbinom(ii - 1, size = popsize, prob = probvec)
}

# Note: we ignore the first column and let almost all elements equal zero, 
# because the probability of sampling k alleles given that you currently have zero
# is zero everywhere except the special case of P00, which is 1 
# (if you have zero now, you will have zero in the next gen with probability 1)
# To address this, we manually set element P00 equal to 1:
tmat[1,1] <- 1

#Now, all the columns should sum to 1:
colSums(tmat)

# Pt <- tmat^2*P0 

# Initial vector of population frequencies:
# Like Buri, suppose we do an experiment with 100 populations
# each with equal allele frequencies (p = q = 0.5) initially.
p_init <- rep(0, popsize + 1)

# note, in this indexing:
#   p_init[1] is the number of populations with 0 alleles
#   p_init[2] is the number with 1 allele
#   p_init[popsize + 1] is the number with 2N alleles

# we want the initial population to have 2N/2 = N p alleles (for 50:50 proportion of alleles)
#  so we fill in the N+1 element of p_init vector with the 
#  number of populations

p_init[N+1] <- 100

## ---- echo = T, eval = F------------------------------------------------------
#  #Now, we will introduce drift in our simulation
#  # after 1 generation of random mating in Wright-Fisher (WF) population,
#  #  the distribution of populations is:
#  tmat %*% p_init
#  
#  # and after 2 generations, it is:
#  tmat %*% (tmat %*% p_init)
#  
#  # and after 3 generations, it is:
#  tmat %*% (tmat %*% (tmat %*% p_init))

## -----------------------------------------------------------------------------
# This is easy enough to iterate.
# Here, we will make a matrix to hold 
# the results, for up to ngen generations

ngen    <- 20
freqmat <- matrix(NA, nrow=popsize+1, ncol=ngen)

curr_pop <- p_init 

for (ii in 1:ngen){
	
	# recompute the new pop frequency distribution
	curr_pop <- tmat %*% curr_pop
	# store it:
	freqmat[,ii] <- curr_pop
}



plot3D::persp3D(x = c(0, probvec), z=freqmat, theta=45, phi=20, contour =F,
                xlab="Allele frequency", zlab="Number of populations",
                ylab="Generations")
 

# Plotting selected generations:
########## 
plot.new()
par(oma=c(1,1,1,1), mfrow=c(2,2), mar=c(2,2,2,0))

pfx <- function(gen){
	plot.new()
	plot.window(xlim=c(0,33), ylim=c(0, 15))
	axis(1, at=seq(-4, 36, by=4))
	axis(2, at=seq(-2, 16, by=2), las=1)
	lines(x=c(16,16), y=c(0, 18), lwd=5, col="gray60")
	title(main = paste("After", gen, ifelse(gen == 1, "gen", "gens")))
}

# Expectation after 1 generation of drift:
pfx(1)
points(x=0:32, y=freqmat[,1], pch=21, bg="red", cex=1.5)

# Expectation after 5 generations of drift
pfx(5)
points(x=0:32, y=freqmat[,5], pch=21, bg="red", cex=1.5)

# Expectation after 10 generations of drift
pfx(10)
points(x=0:32, y=freqmat[,10], pch=21, bg="red", cex=1.5)

# Expectation after 15 generations of drift
pfx(15)
points(x=0:32, y=freqmat[,15], pch=21, bg="red", cex=1.5)

# Restoring old par() configs:
par(oldpar)