#Compute the exact probability of fixation for a new mutation with selection coefficient s
#using the Wright-Fisher model with selection and diploid population size N
fixmut <- function(s = 0.01, N = 10,...) {
#compute and store all post-selection frequencies
p.star <- rep(0,2*N)
for(j in 1:(2*N-1)){
p.star[j] <- j*(1+s)/(j*s + 2*N)
}
#generate entries of the matrix Q
#use the R function that computes the binomial probability density
Q <- matrix(0, nrow = (2*N-1), ncol = (2*N-1))
for(i in 1:(2*N-1)){
for(j in 1:(2*N-1)) {
Q[i, j] <- dbinom(j, 2*N, p.star[i])
}
}
#generate entries fo the vector r (probs. of fixation in one step)
r <- rep(0,(2*N-1))
for(i in 1:(2*N-1)) {
r[i] <- dbinom(2*N, 2*N, p.star[i])
}
#compute u, vector of ultimate fixation probabilities
u <- rep(0,(2*N-1))
I <- diag(2*N-1) #the identity matrix
u <- solve(I - Q, r)
u[1] #output the fixation probability of a new, unique mutation
}
################################
############## 2N = 4 ##########
#neutral mutation
neut <- fixmut(s = 0, N = 2) #should be 1/(2*N)
neut
#advantageous mutation with s = 0.01
fixmut(s = 0.01, N = 2)
fixmut(s = 0.01, N = 2)/neut #fixation prob relative to that of a neutral mutation
#deleterious mutation with s = -0.01
fixmut(s = -0.01, N = 2)
fixmut(s = -0.01, N = 2)/neut #fixation prob relative to that of a neutral mutation
############## 2N = 40 ##########
#neutral mutation
neut <- fixmut(s = 0, N = 20) #should be 1/(2*N)
neut
#advantageous mutation with s = 0.01
fixmut(s = 0.01, N = 20)
fixmut(s = 0.01, N = 20)/neut #fixation prob relative to that of a neutral mutation
#deleterious mutation with s = -0.01
fixmut(s = -0.01, N = 20)
fixmut(s = -0.01, N = 20)/neut #fixation prob relative to that of a neutral mutation
############## 2N = 400 ##########
#neutral mutation
neut <- fixmut(s = 0, N = 200) #should be 1/(2*N)
neut
#advantageous mutation with s = 0.01
fixmut(s = 0.01, N = 200)
fixmut(s = 0.01, N = 200)/neut #fixation prob relative to that of a neutral mutation
#deleterious mutation with s = -0.01
fixmut(s = -0.01, N = 200)
fixmut(s = -0.01, N = 200)/neut #fixation prob relative to that of a neutral mutation