/
drift_variance.R
39 lines (29 loc) · 1.66 KB
/
drift_variance.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
# Visualizing variance due to drift and sampling
total_gens = as.integer(2*65e3/30)
gens = (1:total_gens)
xval = gens*30/1e3
Ne = 10e3
plot( 0, 0 , type="n" , lwd=3 , xlim=c(0,2*65) , ylim=c(0.02,0.08),las=1,bty="n",xlab="Years (thousands)",ylab="Allele Frequency",main="Variance due to neutral drift and sampling")
for( p0 in c(0.05) ) {
samp_0 = 10
samp_t = 100
drift_var = p0 * (1-p0) * ( 1/(2*samp_0) + 1/(2*samp_t) + (gens / (2*Ne))*(1-1/(2*samp_t)) )
polygon( c(xval,rev(xval)) , c(p0 - (drift_var),rev(p0 + (drift_var))) , col=NA , border="#08519c" , lty=2 )
samp_0 = 100
samp_t = 100
drift_var = p0 * (1-p0) * ( 1/(2*samp_0) + 1/(2*samp_t) + (gens / (2*Ne))*(1-1/(2*samp_t)) )
polygon( c(xval,rev(xval)) , c(p0 - (drift_var),rev(p0 + (drift_var))) , col=NA , border="#08519c" , lty=3 )
samp_0 = 10
samp_t = 10
drift_var = p0 * (1-p0) * ( 1/(2*samp_0) + 1/(2*samp_t) + (gens / (2*Ne))*(1-1/(2*samp_t)) )
polygon( c(xval,rev(xval)) , c(p0 - (drift_var),rev(p0 + (drift_var))) , col=NA , border="#08519c" )
drift_var = p0 * (1-p0) * ( gens / (2*Ne) )
polygon( c(xval,rev(xval)) , c(p0 - (drift_var),rev(p0 + (drift_var))) , col="#08519c50" , border=NA )
lines( c(0,tail(xval,n=1)) , c(p0,p0) ,col="#08519c" , lwd=2 , lty=1 )
}
legend("topleft",legend=c("drift alone",
expression(paste('drift, n'[0], ' = ', 10,', n'[t], ' = ', 10)),
expression(paste('drift, n'[0], ' = ', 10,', n'[t], ' = ', 100)),
expression(paste('drift, n'[0], ' = ', 100,', n'[t], ' = ', 100))),
bty="n",
lty=c(1,1,2,3),lwd=c(3,1,1,1),col=c("#08519c50","#08519c","#08519c","#08519c"))