Simulation of population under 20% selection
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
## 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) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment