Simulation of population under 20% selection

20pct_selection.r

## Hardy Weinberg script (3) 

## Simulating changes in genotype frequencies over time 
## under 20% selection against homozygous recessive

## Instruction: Press "-->Source"

##################### CODE: Ignore all of this below ###########################

for (i in 1) {
  parental.alleles <- c(rep("b", times = 100), rep("B", times = 100))
  numbers <- c(sprintf('%0.3d', 1:100), sprintf('%0.3d',1:100))
  x <- paste(parental.alleles, numbers, sep = "")
  BB.f <- NA
  Bb.f <- NA
  bb.f <- NA
  selected.1 <- NA
  k <- 1
  morgue <- NA
  
  while(length(x) > 0) {
    s <- sample(x, replace = FALSE, size = 2)
    a <- paste(s, collapse = "")
    selected.1[k] <- paste(c(substring(a, 1, 1), (substring(a, 5, 5))), collapse = "")
    x <- x[!x %in% s]
    k <- k + 1
  }
  
  BB <- length(selected.1[! selected.1 %in% c('bb', 'Bb', "bB")])
  Bb <- length(selected.1[! selected.1 %in% c('bb', 'BB')])
  bb <- length(selected.1[! selected.1 %in% c('BB', 'Bb', "bB")])
  
  total.a <- BB + Bb + bb
  
  morgue[1] <- floor(bb/5) * 2
  
  BB.f[1] <- BB / total.a
  bb.f[1] <- bb / total.a
  Bb.f[1] <- Bb / total.a
  
  ### GENERATION ONE: 75% of bb survive
  gen1.alleles <- c(rep("b", times = bb*2+Bb-morgue[1]), rep("B", times = BB*2+Bb))
  numbers <- c(sprintf('%0.3d', 1:(bb*2+Bb-morgue[1])), sprintf('%0.3d',1:(BB*2+Bb)))
  x <- paste(gen1.alleles, numbers, sep = "")
  selected.2 <- NA
  k <- 1
  
  while(length(x) > 0) {
    s <- sample(x, replace = FALSE, size = 2)
    a <- paste(s, collapse = "")
    selected.2[k] <- paste(c(substring(a, 1, 1), (substring(a, 5, 5))), collapse = "")
    x <- x[!x %in% s]
    k <- k+1
    
  }
  
  BB <- length(selected.2[! selected.2 %in% c('bb', 'Bb', "bB")])
  Bb <- length(selected.2[! selected.2 %in% c('bb', 'BB')])
  bb <- length(selected.2[! selected.2 %in% c('BB', 'Bb', "bB")])
  
  total.a <- BB + Bb + bb
  
  morgue[2] <- floor(bb/5) * 2
  
  BB.f[2] <- BB / total.a
  bb.f[2] <- bb / total.a
  Bb.f[2] <- Bb / total.a
  
  ### GENERATION TWO: 75% of bb survive
  gen2.alleles <- c(rep("b", times = bb*2+Bb-morgue[2]), rep("B", times = BB*2+Bb))
  numbers <- c(sprintf('%0.3d', 1:(bb*2+Bb-morgue[1])), sprintf('%0.3d',1:(BB*2+Bb)))
  x <- paste(gen2.alleles, numbers, sep = "")
  selected.3 <- NA
  k <- 1
  
  while(length(x) > 0) {
    s <- sample(x, replace = FALSE, size = 2)
    a <- paste(s, collapse = "")
    selected.3[k] <- paste(c(substring(a, 1, 1), (substring(a, 5, 5))), collapse = "")
    x <- x[!x %in% s]
    k <- k+1
  }
  
  BB <- length(selected.3[! selected.3 %in% c('bb', 'Bb', "bB")])
  Bb <- length(selected.3[! selected.3 %in% c('bb', 'BB')])
  bb <- length(selected.3[! selected.3 %in% c('BB', 'Bb', "bB")])
  
  total.a <- BB + Bb + bb
  
  morgue[3] <- floor(bb/5) * 2
  
  BB.f[3] <- BB / total.a
  bb.f[3] <- bb / total.a
  Bb.f[3] <- Bb / total.a
  
  ### GENERATION THREE: 75% of bb survive
  gen3.alleles <- c(rep("b", times = bb*2+Bb-morgue[3]), rep("B", times = BB*2+Bb))
  numbers <- c(sprintf('%0.3d', 1:(bb*2+Bb-morgue[1])), sprintf('%0.3d',1:(BB*2+Bb)))
  x <- paste(gen3.alleles, numbers, sep = "")
  selected.4 <- NA
  k <- 1
  
  while(length(x) > 0) {
    s <- sample(x, replace = FALSE, size = 2)
    a <- paste(s, collapse = "")
    selected.4[k] <- paste(c(substring(a, 1, 1), (substring(a, 5, 5))), collapse = "")
    x <- x[!x %in% s]
    k <- k+1
  }
  
  
  BB <- length(selected.4[! selected.4 %in% c('bb', 'Bb', "bB")])
  Bb <- length(selected.4[! selected.4 %in% c('bb', 'BB')])
  bb <- length(selected.4[! selected.4 %in% c('BB', 'Bb', "bB")])
  
  total.a <- BB + Bb + bb
  
  morgue[4] <- floor(bb/5) * 2
  
  BB.f[4] <- BB / total.a
  bb.f[4] <- bb / total.a
  Bb.f[4] <- Bb / total.a
  
  ### GENERATION FOUR: 75% of bb survive
  gen4.alleles <- c(rep("b", times = bb*2+Bb-morgue[4]), rep("B", times = BB*2+Bb))
  numbers <- c(sprintf('%0.3d', 1:(bb*2+Bb-morgue[1])), sprintf('%0.3d',1:(BB*2+Bb)))
  x <- paste(gen4.alleles, numbers, sep = "")
  selected.5 <- NA
  k <- 1
  
  while(length(x) > 0) {
    s <- sample(x, replace = FALSE, size = 2)
    a <- paste(s, collapse = "")
    selected.5[k] <- paste(c(substring(a, 1, 1), (substring(a, 5, 5))), collapse = "")
    x <- x[!x %in% s]
    k <- k+1
  }
  
  BB <- length(selected.5[! selected.5 %in% c('bb', 'Bb', "bB")])
  Bb <- length(selected.5[! selected.5 %in% c('bb', 'BB')])
  bb <- length(selected.5[! selected.5 %in% c('BB', 'Bb', "bB")])
  
  total.a <- BB + Bb + bb
  
  morgue[5] <- floor(bb/5) * 2
  
  BB.f[5] <- BB / total.a
  bb.f[5] <- bb / total.a
  Bb.f[5] <- Bb / total.a
  
  ### GENERATION FIVE: 75% of bb survive
  gen5.alleles <- c(rep("b", times = bb*2+Bb-morgue[5]), rep("B", times = BB*2+Bb))
  numbers <- c(sprintf('%0.3d', 1:(bb*2+Bb-morgue[1])), sprintf('%0.3d',1:(BB*2+Bb)))
  x <- paste(gen5.alleles, numbers, sep = "")
  selected.6 <- NA
  k <- 1
  
  while(length(x) > 0) {
    s <- sample(x, replace = FALSE, size = 2)
    a <- paste(s, collapse = "")
    selected.6[k] <- paste(c(substring(a, 1, 1), (substring(a, 5, 5))), collapse = "")
    x <- x[!x %in% s]
    k <- k+1
  }
  
  BB <- length(selected.6[! selected.6 %in% c('bb', 'Bb', "bB")])
  Bb <- length(selected.6[! selected.6 %in% c('bb', 'BB')])
  bb <- length(selected.6[! selected.6 %in% c('BB', 'Bb', "bB")])
  
  total.a <- BB + Bb + bb
  
  morgue[6] <- floor(bb/5) * 2
  
  BB.f[6] <- BB / total.a
  bb.f[6] <- bb / total.a
  Bb.f[6] <- Bb / total.a
  
  gen <- seq(1:6)
  generation <- c("P", "1", "2", "3", "4", "5")
  plot(BB.f ~ gen, ylim = c(0, 1), type = "b", col = "blue", 
       ylab = "Genotype frequency", xaxt = "n", 
       xlab = "Generation", pch = 16, lwd = 2,
       main = "Genotype frequencies with 20% selection against bb")
  axis(1, at=1:6, labels=generation)
  lines(Bb.f ~ gen, col = "red", type = "b", pch = 16, lwd = 2)
  lines(bb.f ~ gen, col = "green", type = "b", pch = 16, lwd = 2)
  legend(1, 1, legend=c("BB", "Bb", "bb"),
         col=c("blue", "red", "green"), lty=1, lwd = 2, cex=0.8)
}