Popgen Exercise

Popgen exercise R programs

popgen_drift.q

#############################################################
### Michael Nothnagel, michael.nothnagel@uni-koeln.de ###
### Simulating genetic drift ###
#############################################################

N.s = c(10, 100, 1000, 10000)
n.gen = 50
n.rep = 100
# === Simulate genetic drift === #
for (n in N.s) {
 freqs = matrix(NA, ncol=n.gen+1, nrow=n.rep)
 for (i in 1:n.rep) {
 alleles = c(rep(0,n/2), rep(1,n/2))
 freqs[i,1] = sum(alleles==1) / length(alleles)
 for (j in 1:n.gen) {
 alleles = sample (alleles, length(alleles), replace=T)
 freqs[i,j+1] = sum(alleles==1) / length(alleles)
 }
 }
 assign (paste("freqs.n", n, sep=""), freqs)
 }
summary(freqs.n10 [,n.gen+1])
summary(freqs.n100 [,n.gen+1])
summary(freqs.n1000 [,n.gen+1])
summary(freqs.n10000[,n.gen+1])
# === Graph allele frequency changes === #
pdf("drift_plot.pdf", paper="special", height=4*2, width=4*2, onefile=F)
 split.screen(c(2,2))
 screen(1)
 plot(x=0, y=0, type="n", xlim=c(0,n.gen), ylim=c(0,1), xlab="Generation", ylab="Allele frequency", main="N=10")
 lines (c(0,n.gen), rep(0.5, 2), lty=3, col="#AAAAAA")
 freqs = get(paste("freqs.n", 10 , sep=""))
 for (i in 1:n.rep) {
 lines(0:n.gen,freqs[i,], col="#44AAAA")
 }
 screen(2)
 plot(x=0, y=0, type="n", xlim=c(0,n.gen), ylim=c(0,1), xlab="Generation", ylab="Allele frequency", main="N=100")
 lines (c(0,n.gen), rep(0.5, 2), lty=3, col="#AAAAAA")
 freqs = get(paste("freqs.n", 100 , sep=""))
 for (i in 1:n.rep) {
 lines(0:n.gen,freqs[i,], col="#44AAAA")
 }
 screen(3)
 plot(x=0, y=0, type="n", xlim=c(0,n.gen), ylim=c(0,1), xlab="Generation", ylab="Allele frequency", main="N=1000")
 lines (c(0,n.gen), rep(0.5, 2), lty=3, col="#AAAAAA")
 freqs = get(paste("freqs.n", 1000 , sep=""))
 for (i in 1:n.rep) {
 lines(0:n.gen,freqs[i,], col="#44AAAA")
 }
 screen(4)
 plot(x=0, y=0, type="n", xlim=c(0,n.gen), ylim=c(0,1), xlab="Generation", ylab="Allele frequency", main="N=10000")
 lines (c(0,n.gen), rep(0.5, 2), lty=3, col="#AAAAAA")
 freqs = get(paste("freqs.n", 10000, sep=""))
 for (i in 1:n.rep) {
 lines(0:n.gen,freqs[i,], col="#44AAAA")
 }
 close.screen(all.screens=T)
dev.off()
# === Graph allele frequency after 50 generations === #
pdf("drift_hist.pdf", paper="special", height=4*2, width=4*2, onefile=F)
 split.screen(c(2,2))
 screen(1)
 hist(freqs.n10 [,n.gen+1], main="N=10", xlab="Allele frequency", xlim=c(0,1))
 screen(2)
 hist(freqs.n100 [,n.gen+1], main="N=100", xlab="Allele frequency", xlim=c(0,1))
 screen(3)
 hist(freqs.n1000 [,n.gen+1], main="N=1000", xlab="Allele frequency", xlim=c(0,1))
 screen(4)
 hist(freqs.n10000[,n.gen+1], main="N=10000", xlab="Allele frequency", xlim=c(0,1))
 close.screen(all.screens=T)
dev.off()

popgen.selection.q

#############################################################
### Michael Nothnagel, michael.nothnagel@uni-koeln.de ###
### Calculation allele frequency changes due to selection ###
#############################################################

s.s = c(0.001, 0.01, 0.1)
h = 0.5
n.gen = 100
# === Calculate allele frequency changes === #
for (s in s.s) {
 freqs = rep(NA, n.gen+1)
 af = 0.5
 freqs[1] = 1-af
 for (j in 1:n.gen) {
 omega = 1 - 2*af*(1-af)*h*s - (1-af)*(1-af)*s
 f.het = (1-h*s)*2*af*(1-af)/omega
 f.hom = af*af/omega
 af = f.hom + f.het/2
 freqs[j+1] = 1-af
 }
 assign (paste("freqs.s", s, sep=""), freqs)
 }
# === Report allele frequencies after 100 generations === #
for (s in s.s) {
 cat("s=") ; cat(s) ; cat(": ")
 freqs = get(paste("freqs.s", s, sep=""))
 cat(freqs[n.gen+1]) ; cat("\n")
 }
# === Graph allele frequency changes === #
pdf("selection_plot.pdf", paper="special", height=4*2, width=4*2, onefile=F)
 plot(x=0, y=0, type="n", xlim=c(0,n.gen), ylim=c(0,1), xlab="Generation", ylab="Allele frequency")
 lines (c(0,n.gen), rep(0.5, 2), lty=3, col="#AAAAAA")
 for (s in s.s) {
 freqs = get(paste("freqs.s", s, sep=""))
 lines(0:n.gen,freqs, col="#44AAAA")
 }
dev.off()