Chapter 3
TABLES
"Table
3.gone" - No longer in book as this was deemed to be too heinous.
The Em-algorithm was used to find the allele frequencies on page 56
which are then used in equation 3.4. The EM can be handled very
easily using an Excel spreadsheet. If you still want to see that
go to:
EM algorithm
But you can also use "R" (see P.56 function below, where P.56 refers to the fact that the EM is mentioned on page 56).
Table 3.2 – see the Excel spreadsheet Table 3.2
But you can also use "R" (see Tab3.2 below)
FIGURES, TABLES & Calculations on Specific Pages
All of the stuff below here requires "R" in order to get the job done. See http://www.r-project.org/
. It is with good reason that if you "google" on the letter "R"
the first thing you get to is http://www.r-project.org/.
Anyhow, once you have downloaded "R" and have a "Gui" open, you can
copy and paste the plain font (courier) stuff below.
Fig3.1=function (start=5,stop=50)
{
H=0
for(i in start:stop) H[i-start+1] = 2*choose(i,2)*(1/i)^2
plot(start:stop,H,type='b',xlab='Number of Alleles',
ylab='Frequency of Heterozygotes')
}
Fig3.1()
P.56=function()
{
cat('Enter the number of individuals with type-A blood (=179 in Taylor & Prior)', fill=F ); A = scan(n=1,quiet=T)
cat('Enter the number of individuals with type-B blood (=35 in T & P)', fill=F ); B = scan(n=1,quiet=T)
cat('Enter the number of individuals with type-AB blood (=6 in T & P)', fill=F ); AB = scan(n=1,quiet=T)
cat('Enter the number of individuals with type-O blood (=202 in T & P)', fill=F ); O = scan(n=1,quiet=T)
p=0;q=0;r=0
N=A+B+AB+O
r[1]=sqrt(O/N)
p[1]=(1-r[1])/2
q[1]=p[1]
for(i in 2:10){
AA=p[i-1]^2/(p[i-1]^2+2*p[i-1]*r[i-1])*A
AO=A-AA
BB=q[i-1]^2/(q[i-1]^2+2*q[i-1]*r[i-1])*B
BO=B-BB
p[i]=(2*AA+AO+AB)/(2*N)
q[i]=(2*BB+BO+AB)/(2*N)
r[i]=1-p[i]-q[i]
}
plot(1:10,r,type='b',ylim=c(0,1),xlab='Iteration',ylab='Allele Frequencies',pch='r')
lines(1:10,p,type='b',pch='p')
lines(1:10,q,type='b',pch='q')
return(list(N=N,p=p,q=q,r=r))
}
Tab3.2=function (AA=136,AB=66,BB=17)
{
AA=round(AA,0)
AB=round(AB,0)
BB=round(BB,0)
options(warn=-1)
N=AA+AB+BB
Two.N=2*N
A=2*AA+AB
B=2*BB+AB
f.A=(2*AA+AB)/Two.N
f.B=1-f.A
print(c(AA,AB,BB))
print(c(f.A^2*N,2*f.A*f.B*N,f.B^2*N))
p = function(AA,AB,BB) return(exp(lgamma(N+1)+lgamma(A+1)+lgamma(B+1)+log(2)*AB
-lgamma(AA+1)-lgamma(AB+1)-lgamma(BB+1)-lgamma(Two.N+1)))
start.AA=(A-B)/2
stop.AA=A/2
sto.p=0
obs.p=p(AA,AB,BB)
for(i in start.AA:stop.AA){
t.AA=i
t.AB=A-2*t.AA
t.BB=N-t.AA-t.AB
sto.p[i-start.AA+1]=p(t.AA,t.AB,t.BB)
}
top=max(sto.p)*1.1
plot(c(start.AA-.5,stop.AA+0.5),c(0,top),type='n',xlab='Number of AA genotypes',
ylab='Probability')
tot.p=0
for(i in start.AA:stop.AA){
curr=i-start.AA+1
if(sto.p[curr]<=obs.p) {
tot.p=tot.p+sto.p[curr]
rect(i-.5,0,i+.5,sto.p[curr],density=40)
}
else rect(i-.5,0,i+.5,sto.p[curr])
}
lines(c(AA,AA),c(0,1))
return(round(tot.p,6))
}
Fig3.2=function (recomb=0.5)
{
gam=0
geno=matrix(rep(0,9*100),nc=9)
geno[1,]=c(0,0,.5,0,0,0,.5,0,0)
gam[1]=geno[1,1] + 0.5*geno[1,2] + 0.5*geno[1,4] + 0.25*geno[1,5]
gam[2]=0.5*geno[1,4]+0.25*geno[1,5]+geno[1,7]+0.5*geno[1,8]
gam[3]=0.5*geno[1,2]+geno[1,3]+0.25*geno[1,5]+0.5*geno[1,6]
gam[4]=0.25*geno[1,5]+0.5*geno[1,6]+0.5*geno[1,8]+geno[1,9]
for(i in 2:100){
geno[i,1]=gam[1]^2
geno[i,2]=2*gam[1]*gam[3]
geno[i,3]=gam[3]^2
geno[i,4]=2*gam[1]*gam[2]
geno[i,5]=2*gam[2]*gam[3] + 2*gam[1]*gam[4]
geno[i,6]=2*gam[3]*gam[4]
geno[i,7]=gam[2]^2
geno[i,8]=2*gam[2]*gam[4]
geno[i,9]=gam[4]^2
cross=gam[1]*gam[4]-gam[2]*gam[3]
r=cross/sqrt((gam[1]+gam[3])*(gam[1]+gam[2])*(gam[2]+gam[4])*(gam[1]+gam[2]))
gam[1]=gam[1]-cross*recomb
gam[2]=gam[2]+cross*recomb
gam[3]=gam[3]+cross*recomb
gam[4]=gam[4]-cross*recomb
}
par(mfrow=c(3,3),oma=rep(6,4),mar=rep(0,4))
# First panel
plot(0:99,geno[,1],type='l',ylim=c(0,.5),xlab='Generations',ylab='Genotype Frequencies',xlim=c(0,30),axes=F)
box()
axis(2)
axis(1,labels=F)
text(20,.4,'AABB',cex=2)
# Second panel
plot(0:99,geno[,2],type='l',ylim=c(0,.5),xlab='Generations',ylab='Genotype Frequencies',xlim=c(0,30),axes=F)
box()
axis(1,labels=F)
axis(2,labels=F)
axis(3)
text(20,.4,'AABb',cex=2)
# Third panel
plot(0:99,geno[,3],type='l',ylim=c(0,.5),xlab='Generations',ylab='Genotype Frequencies',xlim=c(0,30),axes=F)
box()
axis(1,labels=F)
axis(2,labels=F)
text(20,.4,'AAbb',cex=2)
plot(0:99,geno[,4],type='l',ylim=c(0,.5),xlab='Generations',ylab='Genotype Frequencies',xlim=c(0,30),axes=F)
box()
axis(1,labels=F)
axis(2,labels=F)
text(20,.4,'AaBB',cex=2)
plot(0:99,geno[,5],type='l',ylim=c(0,.5),xlab='Generations',ylab='Genotype Frequencies',xlim=c(0,30),axes=F)
box()
axis(1,labels=F)
axis(2,labels=F)
text(20,.4,'AaBb',cex=2)
plot(0:99,geno[,6],type='l',ylim=c(0,.5),xlab='Generations',ylab='Genotype Frequencies',xlim=c(0,30),axes=F)
box()
axis(4)
axis(1,labels=F)
text(20,.4,'Aabb',cex=2)
plot(0:99,geno[,7],type='l',ylim=c(0,.5),xlab='Generations',ylab='Genotype Frequencies',xlim=c(0,30),axes=F)
box()
axis(2)
axis(1)
text(20,.4,'aaBB',cex=2)
plot(0:99,geno[,8],type='l',ylim=c(0,.5),xlab='Generations',ylab='Genotype Frequencies',xlim=c(0,30),axes=F)
box()
axis(2,labels=F)
text(20,.4,'aaBb',cex=2)
plot(0:99,geno[,9],type='l',ylim=c(0,.5),xlab='Generations',ylab='Genotype Frequencies',xlim=c(0,30),axes=F)
box()
axis(1)
axis(2,labels=F)
text(20,.4,'aabb',cex=2)
mtext('Generations',1,outer=T,line=3,cex=1.25)
mtext('Genotype Frequencies',2,outer=T,line=3,cex=1.25)
}
Examples
Fig3.2()
Fig3.2(0.1) # To draw Figure 3.4
Fig3.3=function (recomb=.5,add=F,Ngen=15)
{
gam=0
geno=matrix(rep(0,9*100),nc=9)
r=-1
geno[1,]=c(0,0,.5,0,0,0,.5,0,0)
gam[1]=geno[1,1] + 0.5*geno[1,2] + 0.5*geno[1,4] + 0.25*geno[1,5]
gam[2]=0.5*geno[1,4]+0.25*geno[1,5]+geno[1,7]+0.5*geno[1,8]
gam[3]=0.5*geno[1,2]+geno[1,3]+0.25*geno[1,5]+0.5*geno[1,6]
gam[4]=0.25*geno[1,5]+0.5*geno[1,6]+0.5*geno[1,8]+geno[1,9]
for(i in 2:100){
geno[i,1]=gam[1]^2
geno[i,2]=2*gam[1]*gam[3]
geno[i,3]=gam[3]^2
geno[i,4]=2*gam[1]*gam[2]
geno[i,5]=2*gam[2]*gam[3] + 2*gam[1]*gam[4]
geno[i,6]=2*gam[3]*gam[4]
geno[i,7]=gam[2]^2
geno[i,8]=2*gam[2]*gam[4]
geno[i,9]=gam[4]^2
cross=gam[1]*gam[4]-gam[2]*gam[3]
gam[1]=gam[1]-cross*recomb
gam[2]=gam[2]+cross*recomb
gam[3]=gam[3]+cross*recomb
gam[4]=gam[4]-cross*recomb
cross=gam[1]*gam[4]-gam[2]*gam[3]
r[i]=cross/sqrt((gam[1]+gam[3])*(gam[1]+gam[2])*(gam[2]+gam[4])*(gam[1]+gam[2]))
}
par(cex=1.25)
if(add==F){
plot(0:Ngen,r[1:(Ngen+1)],xlab='Generation',ylab='Correlation',
type='l',ylim=c(-1,0))
}
else {lines(0:Ngen,r[1:(Ngen+1)])}
}
Examples
Fig3.3()
Fig3.3(0.1) # To draw figure 3.5
Fig3.6=function (p=0.6)
{
geno=matrix(rep(0,3*16),nc=3)
geno[1,1]=p^2
geno[1,2]=2*p*(1-p)
geno[1,3]=(1-p)^2
F=(2*p*(1-p)-geno[1,2])/(2*p*(1-p))
for(i in 2:16){
geno[i,1]=geno[i-1,1]+0.25*geno[i-1,2]
geno[i,2]=0.5*geno[i-1,2]
geno[i,3]=geno[i-1,3]+0.25*geno[i-1,2]
F[i]=(2*p*(1-p)-geno[i,2])/(2*p*(1-p))
}
plot(0:15,geno[,1],type='l',ylim=c(0,1),
xlab='Generations',ylab='Frequencies (solid)& F (dashed)')
lines(0:15,geno[,2])
lines(0:15,geno[,3])
lines(0:15,F,lty=2,lwd=2)
points(rep(0,3),geno[1,],pch=19)
text(10,p+.05,'AA',cex=1.25,pos=4)
text(10,(1-p)+.05,'BB',cex=1.25,pos=4)
text(10,.05,'AB',cex=1.25,pos=4)
}
Examples
Fig3.6()
Fig3.6(0.75) # To draw different starting conditions
Fig3.9=function (w=c(0.7635,1,0.0001),p.start=c(0.0001,0.9999))
{
p=0.0001
t=0:42
w.AA=w[1]
w.AS=w[2]
w.SS=w[3]
p[1]=p.start[1]
for(i in 2:43){
ave.w=p[i-1]^2*w[1]+2*p[i-1]*(1-p[i-1])*w[2]+(1-p[i-1])^2*w[3]
p[i]=(p[i-1]^2*w[1]+p[i-1]*(1-p[i-1])*w[2])/ave.w
}
plot(t,p,type='l',xlim=c(0,40),ylim=c(0,1),xlab='Generation',
ylab='Allele Frequency (A)')
n.starts=length(p.start)
for(j in 2:n.starts){
p[1]=p.start[j]
for(i in 2:43){
ave.w=p[i-1]^2*w[1]+2*p[i-1]*(1-p[i-1])*w[2]+(1-p[i-1])^2*w[3]
p[i]=(p[i-1]^2*w[1]+p[i-1]*(1-p[i-1])*w[2])/ave.w}
lines(t,p)
}
}
Examples
Fig3.9()
Fig3.09(w=c(1,1,0),p.start=c(0.0001,0.5,0.9)) # Fig. 3.14
Fig3.10=function (w=c(0.7635,1,0.0001),x=c(0,1),y=c(-.1,0.6))
{
p=0.0001
w.AA=w[1]
w.AB=w[2]
w.BB=w[3]
p=seq(0.01,0.99,0.01)
ave.w=p^2*w[1]+2*p*(1-p)*w[2]+(1-p)^2*w[3]
delta.p=p*(1-p)*(p*(w.AA-w.AB)+(1-p)*(w.AB-w.BB))/ave.w
plot(p,delta.p,type='l',xlim=x,ylim=y)
abline(0,0)
delta.p=function(p){
ave.w=p^2*w[1]+2*p*(1-p)*w[2]+(1-p)^2*w[3]
return(p*(1-p)*(p*(w.AA-w.AB)+(1-p)*(w.AB-w.BB))/ave.w)
}
if(sign(delta.p(0.00001))!=sign(delta.p(.99999))){
equil=uniroot(delta.p,c(0.00001,.99999))$root
lines(rep(equil,2),c(-1,1),lty=2)
}
}
Examples
Fig3.10()
Fig3.10(x=c(0.6,1),y=c(-0.025,0.05)) # Zoomed-in version from Fig. 3.11
Fig3.10(c(1,1,0)) # Figure 3.15
Fig3.12=function (w=c(0.7635,1,0.0001),x=c(0,1),y=c(0,1))
{
p=0.0001
w.AA=w[1]
w.AB=w[2]
w.BB=w[3]
p=seq(0.01,0.99,0.01)
ave.w = function(p) return(p^2*w[1]+2*p*(1-p)*w[2]+(1-p)^2*w[3])
plot(p,ave.w(p),type='l',xlim=x,ylim=y)
slope.w=function(p){
return(2*(p*w.AA + (1-2*p)*w.AB + (1-p)*w.BB))
}
d.w=function(p) return(2*(p*(w.AA-w.AB)+(1-p)*(w.AB-w.BB)))
equil=(w.BB-w.AB)/(w.BB+w.AA-2*w.AB)
lines(rep(equil,2),c(0,ave.w(equil)),lty=2)
if(equil<0.0001) {print(noquote('Fixation at p = 0'))
return(noquote(' '))}
if(equil>0.9999) {print(noquote('Fixation at p = 1'))
return(noquote(' '))}
dd=w.AA-2*w.AB+w.BB
if(dd>0) noquote(paste('Unstable equilibrium at p = ',equil))
else(noquote(paste('Stable equilibrium at p = ',equil)))
}
Examples
Fig3.12() # stable equilibrium from original figure
Fig3.12(c(1,0,1))
# Unstable equilibrium
Fig3.12(c(1,1,0)) # lethal recessive (Figure 3.16)
Fig3.12(c(0,1,1)) # lethal dominant
Fig3.13=function (w=c(0.7635,1,0.0001),x=c(.65,.95),y=c(0.77,0.82),
tan.int=.01)
{
w.AA=w[1]
w.AB=w[2]
w.BB=w[3]
p=seq(0.01,0.99,0.01)
ave.w = function(p) return(p^2*w[1]+2*p*(1-p)*w[2]+(1-p)^2*w[3])
slope.w=function(p){
return(2*(p*(w.AA-w.AB)+(1-p)*(w.AB-w.BB)))
}
plot(p,ave.w(p),type='l',xlim=x,ylim=y)
tangent=function(p){
b=slope.w(p)
a=ave.w(p)-p*b
abline(a,b)
}
tans=seq(x[1],x[2],tan.int)
n=length(tans)
for(i in 1:n) tangent(tans[i])
}
Examples
Fig3.13()# Version of Fig. 3.13 with more tangents
Fig3.17=function (N=10)
{
sim.line=function(){
p=0
p[1]=0.5
for(i in 2:100){
p[i]=rbinom(1,2*N,p[i-1])/(2*N)
}
return(p)
}
p=sim.line()
plot(1:100,p,type='l',xlab='Generation',ylab='p',ylim=c(0,1),main=paste('N = ',N))
for(i in 2:100){
p=sim.line()
lines(1:100,p)
}
}
Examples
Fig3.17()
Fig3.17(100) # Fig. 3.18
Fig3.17(1000)# Fig. 3.19
Fig3.20=function (Nids=5,gens=c(5,10,20,30))
{
par(mfrow=c(3,2))
my.hist=function(label){
plot(0:N.alleles,vec,xlim=c(-.5,N.alleles+.5),ylim=c(0,max(vec)),
type='n',axes=F,xlab='Number of A alleles',ylab='Frequency')
if(label==1) {title(main=paste(label,' Generation'))}
else {title(main=paste(label,' Generations'))}
axis(1)
axis(2)
for(i in 1:states){
rect(i-1.5,0,i-.5,vec[i])
}
}
N.alleles=2*Nids
print(noquote(paste(' There are ',N.alleles,' alleles per deme')))
print(noquote("How many alleles should be A alleles at the dawn of time?"))
N.start=scan(n=1)
states=N.alleles+1
Tr=diag(states)
for(i in 1:states){
p=(i-1)/N.alleles
Tr[i,]=dbinom(0:N.alleles,N.alleles,p)
}
vec=c(rep(0,N.start),1,rep(0,states-N.start-1))
my.hist(0)
vec=as.vector(vec%*%Tr)
my.hist(1)
i.count=1
for(i in 1:gens[4]){
vec=as.vector(vec%*%Tr)
if(i==gens[i.count]){
my.hist(gens[i.count])
i.count=i.count+1
}
}
}
Examples
Fig3.20()# And answer 5 of 10 alleles
Fig3.20(20)# And answer 20 of 40 alleles to get Fig. 3.21
Fig3.20(40)# And answer 40 of 80 alleles to get Fig. 3.22
Fig3.20(24)# And answer 36 of 48 alleles to get Fig. 3.23
#WARNING: All of the above examples have multiple panes in the plotting area. To return to one pane:
par(mfrow=c(1,1))
Fig3.24=function (N=5,add=F)
{
Fst=0
for(i in 2:110) Fst[i]=(1/(2*N)+(1-1/(2*N))*Fst[i-1])
if(add==F) {plot(1:110,Fst,ylim=c(0,1),xlim=c(0,100),type='l',xlab='Generations')}
else {lines(1:110,Fst)}
}
Examples
Fig3.24() # First line at N=5
Fig3.24(10,T) # Add second line at N=10
Fig3.24(20,T) # Add second line at N=20
Fig3.24(40,T) # Add second line at N=40
Fig3.25=function (m=0.1,add=F)
{
p1=1
p2=0
for(i in 2:60) {
p1[i]=p1[i-1]*(1-m)+p2[i-1]*m
p2[i]=p2[i-1]*(1-m)+p1[i-1]*m
}
if(add==F){
plot(1:60,p1,ylim=c(0,1),xlim=c(0,50),type='l',
xlab='Generations',ylab='p')
}
else{(lines(1:60,p1))}
lines(1:60,p2)
}
Examples
Fig3.25() # Lines at m = 0.1
Fig3.25(0.05, T) # Add lines at m = 0.05
Fig3.26=function (N=100,n.m=2,add=F)
{
m=n.m/N
Fst=0
Fst.equil=1/(1+4*N*m)
print(noquote(paste('Drift/Migration equilibrium at Fst = ',Fst.equil)))
for(i in 2:210) Fst[i]=(1/(2*N)+(1-1/(2*N))*Fst[i-1])*(1-m)^2
if(add==F){
plot(1:210,Fst,ylim=c(0,1.2*(1/(4*N*m))),xlim=c(0,200),
type='l',xlab='Generations')}
else {
lines(1:210,Fst)}
}
Examples
Fig3.26()
Fig3.26(n.m=1,add=T) # Less migration, so higher equilibrium Fst