#################################### # Chap 3: Fundamentals of genetics # #################################### # Function from previous chapters ################################# meiosis <- function(GF,GM) { from.GM <- rbinom(length(GF),1,0.5) GS <- GF GS[from.GM==1] <- GM[from.GM==1] return(GS) } cross <- function(fa,mo,model) { pat <- meiosis(fa$pat,fa$mat) mat <- meiosis(mo$pat,mo$mat) x <- (pat==model$allele) + (mat==model$allele) m <- model$mu + x*model$alpha + (x==1)*model$delta x <- m + model$sigma*rnorm(length(x)) return(list(pat=pat,mat=mat,pheno=x)) } plot.cross <- function(cross.name, cross) { d <- density(cross$pheno) plot(d,main="",lwd=1.5, xlab = paste("N =",length(cross$pheno),", Bandwidth =",round(d$bw,4))) title(paste(cross.name,": mean = ",round(mean(cross$pheno),2), ", sd = ",round(sd(cross$pheno),2),sep=""),cex.main=1.1) abline(h=0,col=gray(0.3),lwd=1.5) } # End: Function from previous chapters ###################################### N <- 10^5 model <- list(mu=5,alpha=1,delta=0,sigma=0.5,allele="A") pheno1 <- rnorm(N,model$mu,model$sigma) pheno2 <- rnorm(N,model$mu + 2*model$alpha,model$sigma) a <- rep("a",N) A <- rep("A",N) IB1 <- list(pat=a,mat=a,pheno=pheno1) IB2 <- list(pat=A,mat=A,pheno=pheno2) F1 <- cross(IB1,IB2,model) BC1 <- cross(IB1,F1,model) BC2 <- cross(IB2,F1,model) F2 <- cross(F1,F1,model) Fn.fa <- cross(F1,F1,model) Fn.mo <- cross(F1,F1,model) for (g in 2:20) { New.fa <- cross(Fn.fa,Fn.mo,model) New.mo <- cross(Fn.fa,Fn.mo,model) Fn.fa <- New.fa; Fn.mo <- New.mo } RI <- list(pat=c(New.fa$pat,New.mo$pat), mat=c(New.fa$mat,New.mo$mat), pheno=c(New.fa$pheno,New.mo$pheno)) op <- par(mfrow=c(2,3)) plot.cross("IB1",IB1) plot.cross("IB2",IB2) plot.cross("F1",F1) plot.cross("BC2",BC2) plot.cross("F2",F2) plot.cross("RI",RI) par(op) x <- (RI$pat=="A") + (RI$mat=="A") table(x) table(x)/sum(table(x)) Fn.fa <- cross(F1,F1,model) Fn.mo <- cross(F1,F1,model) hetero.average <- mean(c(Fn.fa$pat != Fn.fa$mat, Fn.mo$pat != Fn.mo$mat)) for (g in 2:20) { New.fa <- cross(Fn.fa,Fn.mo,model) New.mo <- cross(Fn.fa,Fn.mo,model) Fn.fa <- New.fa; Fn.mo <- New.mo hetero.average <- c(hetero.average, mean(c(Fn.fa$pat != Fn.fa$mat, Fn.mo$pat != Fn.mo$mat))) } plot(hetero.average,type="l",xlab="generations",lwd=1.5) H <- 0.5^(1:20) lines(H,col=gray(0.4),lwd=1.5) hetero.states <- paste(c(1,2,1,2),"-",c(0,0,1,1),sep="") hetero.states Q <- matrix(0,4,4) Q rownames(Q) <- colnames(Q) <- hetero.states # Transition probabilities Q["1-0",c("1-0","1-1")] <- c(0.5,0.25) Q["2-0","1-1"] <- 1 Q["1-1",] <- c(0.25,0.125,0.25,0.25) Q["2-1",c("1-1","2-1")] <- c(0.25,0.5) Q initial <- c(0,0,1,0) het.count <- c(0.5,0,1,0.5) hetero.prob <- NULL QQ <- Q for (g in 1:20) { hetero.prob[g] <- initial %*% QQ %*% het.count QQ <- QQ %*% Q } lines(hetero.prob,lty=2,lwd=1.5) legend(14.5,0.5, legend=c("simulation","theoretical","selfing"), lty=c(1,2,1),col=gray(c(0,0,0.4)),lwd=rep(1.5,3)) meiosis.rec <- function(GF,GM,rec.frac) { N <- nrow(GF) GS <- GF from.GM <- rbinom(N,1,0.5) GS[from.GM==1,1] <- GM[from.GM==1,1] rec <- rbinom(N,1,rec.frac) from.GM <- from.GM*(1-rec) + (1-from.GM)*rec GS[from.GM==1,2] <- GM[from.GM==1,2] return(GS) } cross.rec <- function(fa, mo,QTL.model,rec.frac) { pat <- meiosis.rec(fa$pat,fa$mat,rec.frac) mat <- meiosis.rec(mo$pat,mo$mat,rec.frac) x <- (pat[,1]==QTL.model$allele) + (mat[,1]==QTL.model$allele) m <- QTL.model$mu + x*QTL.model$alpha + (x==1)*QTL.model$delta y <- m + QTL.model$sigma*rnorm(length(x)) return(list(pat=pat,mat=mat,pheno=y)) } rec.count <- function(Fn) { loc1 <- c(Fn$pat[,1],Fn$mat[,1]) loc2 <- c(Fn$pat[,2],Fn$mat[,2]) cross.tab <- table(loc1,loc2) theta <- (cross.tab[2,1]+cross.tab[1,2])/sum(cross.tab) return(theta) } N <- 10^5 a <- rep("a",N) b <- rep("b",N) A <- rep("A",N) B <- rep("B",N) IB1 <- list(pat=cbind(a,b),mat=cbind(a,b),pheno=pheno1) IB2 <- list(pat=cbind(A,B),mat=cbind(A,B),pheno=pheno2) F1 <- cross.rec(IB1,IB2,model,0) theta.BC <- theta.F2 <- theta.RI <- NULL rec.frac <- seq(0,0.5,by=0.05) for (theta in rec.frac) { BC2 <- cross.rec(IB2,F1,model,theta) theta.BC <- c(theta.BC,rec.count(BC2)) F2 <- cross.rec(F1,F1,model,theta) theta.F2 <- c(theta.F2,rec.count(F2)) Fn.fa <- cross.rec(F1,F1,model,theta) Fn.mo <- cross.rec(F1,F1,model,theta) for (g in 2:20) { New.fa <- cross.rec(Fn.fa,Fn.mo,model,theta) New.mo <- cross.rec(Fn.fa,Fn.mo,model,theta) Fn.fa <- New.fa; Fn.mo <- New.mo } RI <- list(pat=rbind(New.fa$pat,New.mo$pat), mat=rbind(New.fa$mat,New.mo$mat), pheno=c(New.fa$pheno,New.mo$pheno)) theta.RI <- c(theta.RI, rec.count(RI)) } plot(rec.frac,theta.RI,type="l",xlab="theta", ylab="rec. fraction",lwd=1.5) lines(rec.frac,theta.F2,lty=2,lwd=1.5) lines(rec.frac,theta.BC,lty=3,lwd=1.5) legend(0,0.5,legend=c("RI","F2","BC"),lty=1:3,lwd=rep(1.5,3)) lines(rec.frac,4*rec.frac/(1+6*rec.frac),col=gray(0.25),lwd=1.5)