http://statgen.us/index.php?action=history&feed=atom&title=Popgen_Exercise 2026-04-05T20:08:46Z Revision history for this page on the wiki MediaWiki 1.26.2 http://statgen.us/index.php?title=Popgen_Exercise&diff=264&oldid=prev 2017-02-23T18:50:40Z

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

Serveradmin