I just added another function to my evolutionary genetics teaching R package, PopGen to simulate the simultaneous actions of migration, natural selection, & genetic drift.
Here is what the function code looks like:
msd<-function(p0=c(0.5,0.5),Ne=c(100,100),
w=list(c(1,1,1),c(1,1,1)),m=c(0.01,0.01),ngen=400,
colors=c("red","blue"),...){
if(hasArg(show.legend)) show.legend=list(...)$show.legend
else show.legend<-TRUE
w<-lapply(w,function(w) (w/max(w))[3:1])
gametes<-lapply(Ne,function(Ne) rep(0,2*Ne))
gametes<-mapply(function(p0,g,N){
g[1:round(p0*2*N)]<-1
g},p0=p0,g=gametes,N=Ne,SIMPLIFY=FALSE)
p<-lapply(gametes,mean)
for(i in 1:ngen){
genotypes<-lapply(gametes,function(x) matrix(sample(x),
length(x)/2,2))
migrants<-mapply(function(N,m) which(runif(N)<=m),N=Ne,
m=m,SIMPLIFY=FALSE)
for(j in 1:length(genotypes)){
to<-if(j==1) 2 else 1
genotypes[[to]]<-rbind(genotypes[[to]],
genotypes[[j]][migrants[[j]],])
}
for(j in 1:length(genotypes)){
if(length(migrants[[j]])>0)
genotypes[[j]]<-genotypes[[j]][-migrants[[j]],]
}
fitness<-mapply(function(x,w) w[rowSums(x)+1],x=genotypes,
w=w,SIMPLIFY=FALSE)
selected<-mapply(function(prob,N,Ne)
cbind(sample(N,Ne,prob=prob,replace=TRUE),
sample(N,Ne,prob=prob,replace=TRUE)),prob=fitness,
N=sapply(genotypes,nrow),Ne=Ne,SIMPLIFY=FALSE)
copy<-lapply(Ne,function(Ne)
matrix(sample(1:2,2*Ne,replace=TRUE),Ne,2))
gametes<-mapply(function(g,s,c) c(diag(g[s[,1],][,c[,1]]),
diag(g[s[,2],][,c[,2]])),
g=genotypes,s=selected,c=copy,SIMPLIFY=FALSE)
for(j in 1:2) p[[j]][i+1]<-mean(gametes[[j]])
}
plot(0:ngen,p[[1]],type="l",col=colors[1],lwd=2,ylim=c(0,1),
xlab="time (generations)",ylab="f(A)")
lines(x=0:ngen,y=p[[2]],col=colors[2],lwd=2)
if(show.legend) legend(x="topright",legend=1:2,lty=1,col=colors,
lwd=2,bg=make.transparent("white",0.8))
invisible(p)
}
in which the function name, msd
, stands for migration,
selection, & drift. It is a relatively simple modification of
another function, drift.selection
, but in which individuals
are swapped from one population to another at some rate. It's pretty flexible
in that, for instance, we simulate things like asymmetric migration rates
between populations (the function argument m
should contain
emigration rates), differing effective population sizes, and
different genotypic selection between populations.
Here are a few of examples of how it might be used.
First, load package:
library(PopGen)
packageVersion("PopGen")
## [1] '0.9.1'
Now, genetic drift absent migration between populations:
msd(m=c(0,0),Ne=c(200,200),ngen=500) ## with NO migration
msd(m=c(0,0),Ne=c(200,200),ngen=500)
Next, migration but no selection:
Ne<-c(200,200)
msd(m=1/Ne,Ne=Ne,ngen=500) ## relatively low migration
msd(m=10/Ne,Ne=Ne,ngen=500) ## relatively high migration
Now we're ready to introduce divergent natural selection between populations.
First, for reference, absent migration:
msd(p0=c(0.25,0.75),Ne=Ne,w=list(c(1,0.95,0.9),c(0.9,0.95,1)),
m=c(0,0),ngen=500)
Next, with some or a lot of migration:
msd(p0=c(0.25,0.75),Ne=Ne,w=list(c(1,0.95,0.9),c(0.9,0.95,1)),
m=1/Ne,ngen=500) ## some migration
msd(p0=c(0.25,0.75),Ne=Ne,w=list(c(1,0.95,0.9),c(0.9,0.95,1)),
m=10/Ne,ngen=500) ## lots of migration
Finally, what happens when one population is much larger than the other?
msd(p0=c(0.25,0.75),Ne=c(100,1000),w=list(c(1,0.95,0.9),c(0.9,0.95,1)),
m=1/500,ngen=500,show.legend=F)
legend(x="topright",legend=c(expression(paste(N[e],"= 100")),
expression(paste(N[e],"= 1,000"))),lty=1,col=c("red","blue"),
lwd=2,bg=phytools::make.transparent("white", 0.8))
Cool.
Thanks for sharing, nice post! Post really provice useful information!
ReplyDeleteGiaonhan247 chuyên dịch vụ váºn chuyển hà ng Ä‘i Campuchia, váºn chuyển hà ng Ä‘i Là o cÅ©ng như giải đáp gá»i hà ng Ä‘i mỹ bao nhiêu 1kg và chia sẻ bảng giá cước gá»i hà ng Ä‘i mỹ và bảng giá cước phà gá»i hà ng Ä‘i canada, cước phà gá»i hà ng Ä‘i úc giá rẻ nhất.