Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Given a population of AA, Aa, and aa genotypes, how do the two Alleles (A and a) propagate across generations?
rm(list=ls())
library(ggplot2)
library(reshape2)
#Function Definition
#Bookkeeping
calcAllelesInGeneration <- function(x) {
AA = sum(x==1)
AB = sum(x==2)
BB = sum(x==3)
#print(unlist(x))
#print(c(AA, AB, BB))
num.A <- (2 * AA) + AB
num.B <- (2 * BB) + AB
return(c(num.A, num.B))
}
#Given a parent individual, get one of their Alleles in the gamete
getGamete <- function(indiv, prob.big.A) {
if (indiv == 1) return(1) #AA
if (indiv == 3) return(0) #aa
#if Parent is Aa, the gamete is one binomial trial with prob.big.A
if (indiv == 2) return(rbinom(1, size=1, prob.big.A)) #Aa
}
#Two parental Gametes combine to form a zygote
combineGametes <- function(mg,dg) {
if ((mg == 1) && (dg == 1)) return(1) #AA
else if ((mg == 0) && (dg == 0)) return(3) #aa
else return(2) #Aa
}
#Given two individuals, get an offspring for next generation
getOffspring <- function(mom, dad,p) {
momGam <- getGamete(mom,p)
dadGam <- getGamete(dad,p)
return(combineGametes(momGam, dadGam))
}
getNextGen <- function(x, prob) {
nextgen <- list()
for(i in seq(1, kStartPop, by=2)) {
firstborn <- getOffspring(x[i], x[i+1], prob)
secondborn <- getOffspring(x[i], x[i+1], prob)
nextgen[i] <- firstborn
nextgen[i+1] <- secondborn
}
return(unlist(nextgen))
}
simulationOneTrial <- function(start.pop, prob, knumGenerations, trial.index) {
df.allele <- NULL
df.gen <- NULL
#Keep track of the individuals in each generation
df.gen <- rbind(df.gen, start.pop)
x <- start.pop
for ( gen in 1:knumGenerations) {
#print(unlist(nex))
#print(paste("Gen:", gen))
numA <- calcAllelesInGeneration(x)
df.gen <- rbind(df.gen, x)
df.allele <- rbind(df.allele, c(gen, trial.index, numA))
nex <- getNextGen(x, prob)
x <- sample(nex) #shuffle the population for breeding
}
return(as.data.frame(df.allele))
}
###plotting function
plotAllelesWithTime <- function(df, num.trials) {
colorRange<-colorRampPalette(c(rgb(0,0,1), rgb(1,0.7,0) ))
p <- ggplot(df, aes(x= Generation, y= value, group=Trial, color=factor(Trial))) + geom_line()
p <- p + scale_colour_manual(values = colorRange(num.trials),
name="Trial")
p <- p + labs(title = "Allele A Frequencies Across Generations")
p <- p + ylab("Number of \"A\" Allele in the Population")
return(p)
}
#-----------------------------
#RunTime Parameters
kStartPop <- 30
knumGenerations <- 100
knumTrials <- 4
prob.big.A <- 0.45
#-----------------------------
#Uniformly distribute AA, Aa and aa
start.pop <- sample(1:3, kStartPop, replace=TRUE)
df <- NULL
df1<- NULL
for (trial.index in 1: knumTrials) {
df1 <- simulationOneTrial(start.pop, prob.big.A, knumGenerations, trial.index)
df <- rbind(df1, df)
}
names(df) <- c("Generation", "Trial", "A", "a")
#Melt and Plot
mdf <- melt(df, id=c("Generation", "Trial"), measure.vars = "A") #Keep only A allele in the molten data fram
plotAllelesWithTime(mdf, knumTrials)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment