Nothing
## function to compute multiallelic Hardy-Weinberg frequencies
## written by Liam J. Revell 2018
hardy.weinberg<-function(p=c(0.5,0.5),alleles=c("A","a"),as.matrix=FALSE){
if(sum(p)!=1) p<-p/sum(p)
nalleles<-max(length(p),length(alleles))
if(nalleles!=length(p)) p<-rep(1/nalleles,nalleles)
if(nalleles!=length(alleles)) {
if(length(p)<=26) alleles<-LETTERS[1:length(p)]
else alleles<-paste("A",1:length(p),sep="")
}
p<-setNames(p,alleles)
HW<-as.matrix(p)%*%t(as.matrix(p))
if(as.matrix) return(HW)
else {
hw<-vector()
k<-1
for(i in 1:length(alleles)){
for(j in i:length(alleles)){
hw[k]<-if(i==j) HW[i,j] else 2*HW[i,j]
names(hw)[k]<-paste(alleles[c(i,j)],collapse="")
k<-k+1
}
}
return(hw)
}
}
## compute multilocus Hardy-Weinberg in which each individual locus is biallelic
## written by Liam J. Revell 2018
multilocus.hw<-function(nloci=2,p=NULL){
if(is.null(p)) p<-rep(0.5,nloci)
if(length(p)!=nloci) nloci<-length(p)
genotypes<-t(apply(cbind(p,1-p),1,hardy.weinberg))
COMBN<-permutations(n=3,r=nloci,set=T,repeats.allowed=T)
ALLELESp<-LETTERS[1:nloci]
ALLELESq<-letters[1:nloci]
GENOTYPE<-vector()
FREQ<-rep(1,nrow(COMBN))
for(i in 1:nrow(COMBN)){
for(j in 1:nloci){
FREQ[i]<-FREQ[i]*genotypes[j,COMBN[i,j]]
gtype<-if(COMBN[i,j]==1) paste(rep(ALLELESp[j],2),collapse="")
else if(COMBN[i,j]==2) paste(c(ALLELESp[j],ALLELESq[j]),collapse="")
else if(COMBN[i,j]==3) paste(rep(ALLELESq[j],2),collapse="")
GENOTYPE[i]<-if(j==1) gtype else paste(GENOTYPE[i],gtype,sep="")
}
}
hw<-setNames(FREQ,GENOTYPE)
hw
}
## function to compute relative frequencies of a phenotypic trait
## for a polygenic trait under a simple additive genetic model
## written by Liam J. Revell 2018, 2019
phenotype.freq<-function(nloci=6,p=NULL,effect=1/nloci){
if(is.null(p)) p<-rep(0.5,nloci)
if(length(effect)==1) effect<-rep(effect,nloci)
else if(length(effect)>1&&length(effect)!=nloci){
cat("The length of \'effect\' should match \'nloci\'.\n")
cat("Setting to recycle first value of \'effect\'.\n")
effect<-rep(effect[1],nloci)
}
genotypes<-t(apply(cbind(p,1-p),1,hardy.weinberg))
COMBN<-permutations(n=3,r=nloci,set=T,repeats.allowed=T)
PHEN<-round(-rowSums((COMBN-2)%*%effect),12)
FREQ<-rep(1,nrow(COMBN))
for(i in 1:nrow(COMBN)){
for(j in 1:nloci) FREQ[i]<-FREQ[i]*genotypes[j,COMBN[i,j]]
}
phen<-sort(unique(PHEN))
if(length(phen)<0.2*length(PHEN)){
type<-"discrete"
freq<-rep(0,length(phen))
for(i in 1:length(phen)) freq[i]<-sum(FREQ[which(PHEN==phen[i])])
plot(phen,freq,type="b",pch=21,bg="grey",
xlab="phenotypic trait value",ylab="relative frequency",
cex=1.5,xlim=range(phen)+c(-0.5,0.5)*mean(effect),
ylim=c(0,max(freq)))
by<-min(rowSums(cbind(-phen[2:length(phen)-1],phen[2:length(phen)])))
for(i in 1:length(phen))
rect(phen[i]-0.4*by,0,phen[i]+0.4*by,
freq[i],border="grey",
col=make.transparent("blue",0.2))
} else {
type<-"continuous"
h<-hist(PHEN,breaks=seq(min(PHEN),max(PHEN),
by=diff(range(PHEN))/20),plot=FALSE)
phen<-h$mids
freq<-h$counts/sum(h$counts)
plot(phen,freq,type="b",pch=21,bg="grey",
xlab="phenotypic trait value",ylab="relative frequency",
cex=1.5,xlim=range(phen)+c(-0.5,0.5)*mean(effect),
ylim=c(0,max(freq)))
by<-phen[2]-phen[1]
for(i in 1:length(phen))
rect(phen[i]-0.5*by,0,phen[i]+0.5*by,
freq[i],border="grey",
col=make.transparent("blue",0.2))
}
object<-list(phenotype=phen,frequency=freq,
nloci=nloci,p=p,effect=effect,
type=type)
class(object)<-"phenotype.freq"
invisible(object)
}
print.phenotype.freq<-function(x,...){
cat("\nObject of class \"phenotype.freq\" containing the frequencies of each value")
cat(paste("\nof a hypothetical polygenic trait determined by the additive effect of",x$nloci))
cat("\ngenetic loci with the following frequencies,\n")
cat(paste(" p(A): ",paste(round(x$p,2),collapse=", "),"\n",sep=""))
cat("and the following additive effect of allelic subsitution,\n")
cat(paste(" a: ",paste(round(x$effect,2),collapse=", "),"\n\n",sep=""))
cat("\nTo plot enter plot(\'object_name\') at the command line interface.\n\n")
}
plot.phenotype.freq<-function(x,...){
phen<-x$phenotype
freq<-x$frequency
if(x$type=="discrete"){
plot(phen,freq,type="b",pch=21,bg="grey",
xlab="phenotypic trait value",ylab="relative frequency",
cex=1.5,xlim=range(phen)+c(-0.5,0.5)*mean(x$effect),
ylim=c(0,max(freq)))
by<-min(rowSums(cbind(-phen[2:length(phen)-1],
phen[2:length(phen)])))
for(i in 1:length(phen))
rect(phen[i]-0.4*by,0,phen[i]+0.4*by,
freq[i],border="grey",
col=make.transparent("blue",0.2))
} else {
plot(phen,freq,type="b",pch=21,bg="grey",
xlab="phenotypic trait value",ylab="relative frequency",
cex=1.5,xlim=range(phen)+c(-0.5,0.5)*mean(x$effect),
ylim=c(0,max(freq)))
by<-phen[2]-phen[1]
for(i in 1:length(phen))
rect(phen[i]-0.5*by,0,phen[i]+0.5*by,
freq[i],border="grey",
col=make.transparent("blue",0.2))
}
}
## function to conduct numerical analysis of phenotypic selection on
## a polygenic quantitative trait
## written by Liam J. Revell 2018
phenotype.selection<-function(nloci=6,p=NULL,effect=1/nloci,beta=0.1,ngen=20,...){
if(hasArg(sleep)) sleep<-list(...)$sleep
else sleep<-0.1
if(is.null(p)) p<-rep(0.5,nloci)
COMBN<-permutations(n=3,r=nloci,set=T,repeats.allowed=T)
PHEN<--rowSums(COMBN-2)*effect
w<-beta*(PHEN-min(PHEN))+1
het<-as.matrix(apply(COMBN,2,function(x) which(x==2)))
homo<-as.matrix(apply(COMBN,2,function(x) which(x==1)))
for(i in 1:ngen){
genotypes<-t(apply(cbind(p,1-p),1,hardy.weinberg))
FREQ<-rep(1,nrow(COMBN))
for(j in 1:nrow(COMBN)){
for(k in 1:nloci) FREQ[j]<-FREQ[j]*genotypes[k,COMBN[j,k]]
}
phen<-unique(PHEN)
freq<-rep(0,length(phen))
for(j in 1:length(phen)) freq[j]<-sum(FREQ[which(PHEN==phen[j])])
dev.hold()
plot(phen,freq,type="b",pch=21,bg="grey",
xlab="phenotypic trait value",ylab="relative frequency",
cex=1.5,xlim=range(phen)+c(-0.5,0.5)*effect,ylim=c(0,max(freq)))
for(j in 1:length(phen))
rect(phen[j]-0.4*effect,0,phen[j]+0.4*effect,
freq[j],border="grey",
col=make.transparent("blue",0.2))
dev.flush()
FREQp<-FREQ*w/sum(FREQ*w)
pp<-vector()
for(j in 1:nloci) pp[j]<-sum(FREQp[het[,j]]/2)+sum(FREQp[homo[,j]])
p<-pp
Sys.sleep(sleep)
}
}
## plot a coalescent genealogy
## written by Liam J. Revell 2018, 2019
coalescent.plot<-function(n=10,ngen=20,colors=NULL,...){
if(hasArg(sleep)) sleep<-list(...)$sleep
else sleep<-0.2
if(hasArg(lwd)) lwd<-list(...)$lwd
else lwd<-2
if(hasArg(mar)) mar<-list(...)$mar
else mar<-c(2.1,4.1,2.1,1.1)
if(is.null(colors)){
colors<-rainbow(n=n)
if(hasArg(col.order)) col.order<-list(...)$col.order
else col.order<-"sequential"
if(col.order=="alternating"){
if(n%%2==1)
ii<-as.vector(rbind(1:((n+1)/2),1:((n+1)/2)+(n+1)/2))
else
ii<-as.vector(rbind(1:(n/2),n/2+1:(n/2)))
colors<-colors[ii]
}
}
popn<-matrix(NA,ngen+1,n)
parent<-matrix(NA,ngen,n)
popn[1,]<-1:n
for(i in 1:ngen){
parent[i,]<-sort(sample(1:n,replace=TRUE))
popn[i+1,]<-popn[i,parent[i,]]
}
plot.new()
par(mar=mar)
plot.window(xlim=c(0.5,n+0.5),ylim=c(ngen,0))
axis(2)
title(ylab="time (generations)")
cx.pt<-2*25/max(n,ngen)
points(1:n,rep(0,n),bg=colors,pch=21,cex=cx.pt)
for(i in 1:ngen){
dev.hold()
for(j in 1:n){
lines(c(parent[i,j],j),c(i-1,i),lwd=lwd,
col=colors[popn[i+1,j]])
}
points(1:n,rep(i-1,n),col="grey",bg=colors[popn[i,]],pch=21,
cex=cx.pt)
points(1:n,rep(i,n),col="grey",bg=colors[popn[i+1,]],pch=21,
cex=cx.pt)
dev.flush()
Sys.sleep(sleep)
}
object<-list(allele=popn,parent=parent)
class(object)<-"coalescent.plot"
invisible(object)
}
print.coalescent.plot<-function(x,...){
cat("\nObject of class \"coalescent.plot\" consisting of a simulated process of allelic\n")
cat(paste("coalescence over",nrow(x$allele)-1,"generations, in a population containing",
ncol(x$allele),"individuals.\n\n"))
cat("The object consists of:\n")
cat(paste(" (1) a ",nrow(x$allele)," x ",ncol(x$allele),
" numeric matrix containing the unique \'alleles\'\n",sep=""))
cat(paste(" present in the population from time=0 to time=",nrow(x$allele)-1,".\n",sep=""))
cat(paste(" (2) a ",nrow(x$parent)," x ",ncol(x$parent),
" numeric matrix giving the parent/offspring\n",sep=""))
cat(paste(" relationships across all ",nrow(x$parent),
" generations of the simulation.\n\n",sep=""))
cat("To re-plot type plot(object_name) at the command line.\n\n")
}
plot.coalescent.plot<-function(x,...){
popn<-x$allele
parent<-x$parent
n<-ncol(popn)
ngen<-nrow(parent)
if(hasArg(sleep)) sleep<-list(...)$sleep
else sleep<-0.2
if(hasArg(lwd)) lwd<-list(...)$lwd
else lwd<-2
if(hasArg(mar)) mar<-list(...)$mar
else mar<-c(2.1,4.1,2.1,1.1)
if(hasArg(colors)) colors<-list(...)$colors
else colors<-NULL
if(is.null(colors)){
colors<-rainbow(n=n)
if(hasArg(col.order)) col.order<-list(...)$col.order
else col.order<-"sequential"
if(col.order=="alternating"){
if(n%%2==1)
ii<-as.vector(rbind(1:((n+1)/2),1:((n+1)/2)+(n+1)/2))
else
ii<-as.vector(rbind(1:(n/2),n/2+1:(n/2)))
colors<-colors[ii]
}
}
plot.new()
par(mar=mar)
plot.window(xlim=c(0.5,n+0.5),ylim=c(ngen,0))
axis(2)
title(ylab="time (generations)")
cx.pt<-2*25/max(n,ngen)
points(1:n,rep(0,n),bg=colors,pch=21,cex=cx.pt)
for(i in 1:ngen){
dev.hold()
for(j in 1:n){
lines(c(parent[i,j],j),c(i-1,i),lwd=lwd,
col=colors[popn[i+1,j]])
}
points(1:n,rep(i-1,n),col="grey",bg=colors[popn[i,]],pch=21,
cex=cx.pt)
points(1:n,rep(i,n),col="grey",bg=colors[popn[i+1,]],pch=21,
cex=cx.pt)
dev.flush()
Sys.sleep(sleep)
}
}
## function to simulate drift & selection
drift.selection<-function(p0=0.5,Ne=100,w=c(1,1,1),ngen=400,nrep=10,
colors=NULL,...){
if(is.null(colors)) colors<-rainbow(nrep)
w<-(w/max(w))[3:1]
gametes<-rep(0,2*Ne)
if(p0>0) gametes[1:round(p0*2*Ne)]<-1
gametes<-replicate(nrep,gametes,simplify=FALSE)
p<-lapply(gametes,mean)
for(i in 1:ngen){
genotypes<-lapply(gametes,function(x) matrix(sample(x),
length(x)/2,2))
fitness<-lapply(genotypes,function(x,w) w[rowSums(x)+1],w=w)
selected<-lapply(fitness,function(prob,x)
cbind(sample(x,prob=prob,replace=TRUE),
sample(x,prob=prob,replace=TRUE)),x=Ne)
copy<-replicate(nrep,matrix(sample(1:2,2*Ne,replace=TRUE),
Ne,2),simplify=FALSE)
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:nrep) 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)")
if(nrep>1)
nulo<-mapply(lines,x=replicate(nrep-1,0:ngen,simplify=FALSE),
y=p[2:nrep],col=colors[2:nrep],lwd=2)
class(p)<-"drift.selection"
attr(Ne,"Ne")<-Ne
attr(p,"w")<-w[3:1]
invisible(p)
}
print.drift.selection<-function(x,...){
ngen<-length(x[[1]])-1
nsim<-length(x)
cat("\nObject of class \"drift.selection\" consisting of")
cat(paste("\nallele frequencies through time from ",nsim," simulation(s) of\n",sep=""))
cat(paste(ngen," generations of genetic drift & natural selection",sep=""))
cat("\nwith normalized fitnesses of:\n")
cat(paste(" W(AA) =",round(attr(x,"w")[1],2),"\n"))
cat(paste(" W(Aa) =",round(attr(x,"w")[2],2),"\n"))
cat(paste(" W(aa) =",round(attr(x,"w")[3],2),"\n"))
cat("\nTo plot enter plot(\'object_name\') at the command line\ninterface.\n\n")
}
plot.drift.selection<-function(x,...){
ngen<-length(x[[1]])-1
nrep<-length(x)
if(hasArg(colors)) colors<-list(...)$colors
else colors<-rainbow(nrep)
if(hasArg(type)) type<-list(...)$type
else type<-"l"
if(hasArg(lwd)) lwd<-list(...)$lwd
else lwd<-2
plot(0:ngen,x[[1]],type=type,col=colors[1],lwd=lwd,ylim=c(0,1),
xlab="time (generations)",ylab="f(A)")
if(nrep>1)
nulo<-mapply(lines,x=replicate(nrep-1,0:ngen,simplify=FALSE),
y=x[2:nrep],col=colors[2:nrep],type=type,lwd=lwd)
}
## function to simulate migration, selection, & drift
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))
attr(p,"p0")<-p0
attr(p,"Ne")<-Ne
attr(p,"w")<-lapply(w,function(w) w[3:1])
attr(p,"m")<-m
class(p)<-"msd"
invisible(p)
}
print.msd<-function(x,...){
cat("\nObject of class \"msd\" containing the results from a numerical simulation")
cat("\nof drift, selection, and migration within and between two populations.")
cat("\n\nThe following parameters were used in the simulation:\n")
cat(paste("\t",paste("Ne[",1:2,"]",sep="",collapse="\t"),"\n",sep=""))
cat(paste("\t",paste(attr(x,"Ne"),collapse="\t")))
cat("\n")
cat("\tm[1->2]\tm[2->1]\n")
cat(paste("\t",paste(attr(x,"m"),collapse="\t")))
cat("\n")
print(matrix(c(attr(x,"w")[[1]],attr(x,"w")[[2]]),2,3,byrow=TRUE,
dimnames=list(c("pop'n 1","pop'n 2"),c("w(AA)","w(Aa)","w(aa)"))))
cat("\n")
cat("\nTo plot enter plot(\'object_name\') at the command line interface.\n\n")
}
plot.msd<-function(x,...){
if(hasArg(show.legend)) show.legend=list(...)$show.legend
else show.legend<-TRUE
if(hasArg(colors)) colors<-list(...)$colors
else colors<-c("red","blue")
if(hasArg(lwd)) lwd<-list(...)$lwd
else lwd<-2
if(hasArg(type)) type<-list(...)$type
else type<-"l"
ngen<-length(x[[1]])-1
plot(0:ngen,x[[1]],type=type,col=colors[1],lwd=lwd,ylim=c(0,1),
xlab="time (generations)",ylab="f(A)")
lines(x=0:ngen,y=x[[2]],col=colors[2],lwd=lwd,type=type)
if(show.legend) legend(x="topright",legend=1:2,lty=1,col=colors,
lwd=lwd,bg=make.transparent("white",0.8))
}
## function to illustrate the central limit theorem (CLT)
clt<-function(nvar=1,nobs=1000,df=c("normal","uniform","exponential","binomial"),
theta=NULL,breaks="Sturges",show=c("sum","mean")){
df<-df[1]
show<-show[1]
if(is.null(theta))
theta<-if(df%in%c("normal","uniform","exponential")) 1 else 0.5
foo<-if(df=="normal") function(n,theta){
rnorm(n,sd=sqrt(theta))
} else if(df=="uniform") function(n,theta){
runif(n,0,theta)
} else if(df=="exponential") function(n,theta){
rexp(n,rate=theta)
} else if(df=="binomial") function(n,theta){
rbinom(n,1,theta)
}
X<-matrix(foo(nvar*nobs,theta),nobs,nvar)
x<-if(show=="sum") rowSums(X) else rowMeans(X)
if(is.integer(breaks)) breaks<-seq(min(x),max(x),by=diff(range(x))/(breaks-1))
obj<-hist(x,border="darkgrey",col=phytools::make.transparent("blue",0.1),
main=paste("CLT: the",show,"of",nvar,df,"distribution(s)"),
breaks=breaks)
lines(obj$mids,obj$counts,type="b",pch=21,bg="grey")
object<-list(dist=obj,data=X,df=df,show=show)
class(object)<-"clt"
invisible(object)
}
## S3 methods for clt
print.clt<-function(x,...){
cat("\nObject of class \"clt\" consisting of:\n")
cat(paste(" (1) ",ncol(x$data)," ",x$df,
"ly distributed independent random variables, each with\n",sep=""))
cat(paste(" (2)",nrow(x$data),"observations (stored in x$data), and\n"))
cat(paste(" (3) a histogram giving the distribution of their observation-wise ",
x$show,".\n\n",sep=""))
}
plot.clt<-function(x,...){
plot(x$dist,border="darkgrey",col=phytools::make.transparent("blue",0.1),
main=paste("CLT: the",x$show,"of",ncol(x$data),x$df,"distribution(s)"))
lines(x$dist$mids,x$dist$counts,type="b",pch=21,bg="grey")
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.