6,689 bytes added,
18:50, 23 February 2017 __NOTITLE__
==Popgen exercise R programs==
popgen_drift.q
#############################################################<br data-attributes="%20/" />### Michael Nothnagel, michael.nothnagel@uni-koeln.de ###<br data-attributes="%20/" />### Simulating genetic drift ###<br data-attributes="%20/" />#############################################################
<br data-attributes="%20/" />N.s = c(10, 100, 1000, 10000)<br data-attributes="%20/" />n.gen = 50<br data-attributes="%20/" />n.rep = 100
# === Simulate genetic drift === #<br data-attributes="%20/" />for (n in N.s) {<br data-attributes="%20/" /> freqs = matrix(NA, ncol=n.gen+1, nrow=n.rep)<br data-attributes="%20/" /> for (i in 1:n.rep) {<br data-attributes="%20/" /> alleles = c(rep(0,n/2), rep(1,n/2))<br data-attributes="%20/" /> freqs[i,1] = sum(alleles==1) / length(alleles)<br data-attributes="%20/" /> for (j in 1:n.gen) {<br data-attributes="%20/" /> alleles = sample (alleles, length(alleles), replace=T)<br data-attributes="%20/" /> freqs[i,j+1] = sum(alleles==1) / length(alleles)<br data-attributes="%20/" /> }<br data-attributes="%20/" /> }<br data-attributes="%20/" /> assign (paste("freqs.n", n, sep=""), freqs)<br data-attributes="%20/" /> }
summary(freqs.n10 [,n.gen+1])<br data-attributes="%20/" />summary(freqs.n100 [,n.gen+1])<br data-attributes="%20/" />summary(freqs.n1000 [,n.gen+1])<br data-attributes="%20/" />summary(freqs.n10000[,n.gen+1])
# === Graph allele frequency changes === #<br data-attributes="%20/" />pdf("drift_plot.pdf", paper="special", height=4*2, width=4*2, onefile=F)<br data-attributes="%20/" /> split.screen(c(2,2))<br data-attributes="%20/" /> screen(1)<br data-attributes="%20/" /> plot(x=0, y=0, type="n", xlim=c(0,n.gen), ylim=c(0,1), xlab="Generation", ylab="Allele frequency", main="N=10")<br data-attributes="%20/" /> lines (c(0,n.gen), rep(0.5, 2), lty=3, col="#AAAAAA")<br data-attributes="%20/" /> freqs = get(paste("freqs.n", 10 , sep=""))<br data-attributes="%20/" /> for (i in 1:n.rep) {<br data-attributes="%20/" /> lines(0:n.gen,freqs[i,], col="#44AAAA")<br data-attributes="%20/" /> }<br data-attributes="%20/" /> screen(2)<br data-attributes="%20/" /> plot(x=0, y=0, type="n", xlim=c(0,n.gen), ylim=c(0,1), xlab="Generation", ylab="Allele frequency", main="N=100")<br data-attributes="%20/" /> lines (c(0,n.gen), rep(0.5, 2), lty=3, col="#AAAAAA")<br data-attributes="%20/" /> freqs = get(paste("freqs.n", 100 , sep=""))<br data-attributes="%20/" /> for (i in 1:n.rep) {<br data-attributes="%20/" /> lines(0:n.gen,freqs[i,], col="#44AAAA")<br data-attributes="%20/" /> }<br data-attributes="%20/" /> screen(3)<br data-attributes="%20/" /> plot(x=0, y=0, type="n", xlim=c(0,n.gen), ylim=c(0,1), xlab="Generation", ylab="Allele frequency", main="N=1000")<br data-attributes="%20/" /> lines (c(0,n.gen), rep(0.5, 2), lty=3, col="#AAAAAA")<br data-attributes="%20/" /> freqs = get(paste("freqs.n", 1000 , sep=""))<br data-attributes="%20/" /> for (i in 1:n.rep) {<br data-attributes="%20/" /> lines(0:n.gen,freqs[i,], col="#44AAAA")<br data-attributes="%20/" /> }<br data-attributes="%20/" /> screen(4)<br data-attributes="%20/" /> plot(x=0, y=0, type="n", xlim=c(0,n.gen), ylim=c(0,1), xlab="Generation", ylab="Allele frequency", main="N=10000")<br data-attributes="%20/" /> lines (c(0,n.gen), rep(0.5, 2), lty=3, col="#AAAAAA")<br data-attributes="%20/" /> freqs = get(paste("freqs.n", 10000, sep=""))<br data-attributes="%20/" /> for (i in 1:n.rep) {<br data-attributes="%20/" /> lines(0:n.gen,freqs[i,], col="#44AAAA")<br data-attributes="%20/" /> }<br data-attributes="%20/" /> close.screen(all.screens=T)<br data-attributes="%20/" />dev.off()
# === Graph allele frequency after 50 generations === #<br data-attributes="%20/" />pdf("drift_hist.pdf", paper="special", height=4*2, width=4*2, onefile=F)<br data-attributes="%20/" /> split.screen(c(2,2))<br data-attributes="%20/" /> screen(1)<br data-attributes="%20/" /> hist(freqs.n10 [,n.gen+1], main="N=10", xlab="Allele frequency", xlim=c(0,1))<br data-attributes="%20/" /> screen(2)<br data-attributes="%20/" /> hist(freqs.n100 [,n.gen+1], main="N=100", xlab="Allele frequency", xlim=c(0,1))<br data-attributes="%20/" /> screen(3)<br data-attributes="%20/" /> hist(freqs.n1000 [,n.gen+1], main="N=1000", xlab="Allele frequency", xlim=c(0,1))<br data-attributes="%20/" /> screen(4)<br data-attributes="%20/" /> hist(freqs.n10000[,n.gen+1], main="N=10000", xlab="Allele frequency", xlim=c(0,1))<br data-attributes="%20/" /> close.screen(all.screens=T)<br data-attributes="%20/" />dev.off()
popgen.selection.q
#############################################################<br data-attributes="%20/" />### Michael Nothnagel, michael.nothnagel@uni-koeln.de ###<br data-attributes="%20/" />### Calculation allele frequency changes due to selection ###<br data-attributes="%20/" />#############################################################
<br data-attributes="%20/" />s.s = c(0.001, 0.01, 0.1)<br data-attributes="%20/" />h = 0.5<br data-attributes="%20/" />n.gen = 100
# === Calculate allele frequency changes === #<br data-attributes="%20/" />for (s in s.s) {<br data-attributes="%20/" /> freqs = rep(NA, n.gen+1)<br data-attributes="%20/" /> af = 0.5<br data-attributes="%20/" /> freqs[1] = 1-af<br data-attributes="%20/" /> for (j in 1:n.gen) {<br data-attributes="%20/" /> omega = 1 - 2*af*(1-af)*h*s - (1-af)*(1-af)*s<br data-attributes="%20/" /> f.het = (1-h*s)*2*af*(1-af)/omega<br data-attributes="%20/" /> f.hom = af*af/omega<br data-attributes="%20/" /> af = f.hom + f.het/2<br data-attributes="%20/" /> freqs[j+1] = 1-af<br data-attributes="%20/" /> }<br data-attributes="%20/" /> assign (paste("freqs.s", s, sep=""), freqs)<br data-attributes="%20/" /> }
# === Report allele frequencies after 100 generations === #<br data-attributes="%20/" />for (s in s.s) {<br data-attributes="%20/" /> cat("s=") ; cat(s) ; cat(": ")<br data-attributes="%20/" /> freqs = get(paste("freqs.s", s, sep=""))<br data-attributes="%20/" /> cat(freqs[n.gen+1]) ; cat("\n")<br data-attributes="%20/" /> }
# === Graph allele frequency changes === #<br data-attributes="%20/" />pdf("selection_plot.pdf", paper="special", height=4*2, width=4*2, onefile=F)<br data-attributes="%20/" /> plot(x=0, y=0, type="n", xlim=c(0,n.gen), ylim=c(0,1), xlab="Generation", ylab="Allele frequency")<br data-attributes="%20/" /> lines (c(0,n.gen), rep(0.5, 2), lty=3, col="#AAAAAA")<br data-attributes="%20/" /> for (s in s.s) {<br data-attributes="%20/" /> freqs = get(paste("freqs.s", s, sep=""))<br data-attributes="%20/" /> lines(0:n.gen,freqs, col="#44AAAA")<br data-attributes="%20/" /> }<br data-attributes="%20/" />dev.off()