fst_dots.R

par(mfrow=c(1,5))
par(oma=c(0.5,0.5,0.5,0.5))
par(mar=c(0.5,0.5,1,0.5))

# fst parameters
# frequency in one population
p_b = 0.2
# scan of frequencies in the second population
p_a = seq(0.01,0.99,by=0.01)
# estimate using Hudson's Fst
fst = ((p_a - p_b)^2) / ( p_a*(1-p_b) + p_b*(1-p_a))
# colors
clr = c("#ffffff","#a6bddb","#1c9099")
  
# plot the different Fst panels
for ( cur_fst in c(0.01,0.10,0.20,0.50,0.75) ) {
  # find the frequency in the second population that produced the closest fst
  cur_p_a = p_a[ which.min( abs(fst - cur_fst) ) ]
  
  # make x and y values and colors
  N = 40
  M = 15
  ypos = c(1:20,22:41)
  mat_x = matrix(NA,nrow=N,ncol=M)
  mat_y = matrix(NA,nrow=N,ncol=M)
  mat_col = matrix(0,nrow=N,ncol=M)
  for ( i in 1:M ) {
    mat_y[,i] = ypos
    mat_x[,i] = i
  }
  
  # sample alleles
  mat_col[1:20,] = rbinom(M*20,2,cur_p_a)
  if ( cur_fst == 0.01 ) {
    # fix alleles in one of the subpopulations
    main_pop = rbinom(M*20,2,p_b)
  }
  mat_col[21:40,] = main_pop
  
  # plot the panel
  plot(t(mat_x),t(mat_y),xlim=c(0,M+1),col=clr[t(mat_col)+1],pch=19,bty="n",yaxt="n",xaxt="n",main=bquote('F'[ST] == .(cur_fst) ),xlab="",ylab="")
  points(t(mat_x),t(mat_y),col=clr[3])
  abline(h=21,lty=3)
}