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
You can’t perform that action at this time.