################################################# # Chap 2: Introduction to experimental genetics # ################################################# n <- 6; p <- 0.5; x <- rbinom(n,2,p) mu <- 5; alpha <- 1; delta <- -1 mu + alpha*x + delta*(x==1) sig <- 0.5 mu + alpha*x + delta*(x==1) + sig*rnorm(n) rnorm(n, mu + alpha*x + delta*(x==1), sig) mu + alpha*x + delta*(x==1) + rnorm(n, sd=sig) par(mfcol=c(1,2)) n <- 10^5; p <- 0.5; x <- rbinom(n,2,p) mu <- 5; alpha <- 1; delta <- -1; sig <- 0.5 y <- mu + alpha*x + delta*(x==1) + sig*rnorm(n) h2 <- var(alpha*x + delta*(x==1))/var(y) d <- density(y) plot(d,main= paste("A recessive model:\n h^2 = ",round(h2,3),sep=""), xlab = paste("N = 100000, Bandwidth =",round(d$bw,4)),lwd=1.5) abline(h=0,col=gray(0.3),lwd=1.5) mu <- 5; alpha <- 1; delta <- -1; sig <- 1 y <- mu + alpha*x + delta*(x==1) + sig*rnorm(n) h2 <- var(alpha*x + delta*(x==1))/var(y) d <- density(y) plot(d,main= paste("A recessive model:\n h^2 = ",round(h2,3),sep=""), xlab = paste("N = 100000, Bandwidth =",round(d$bw,4)), lwd=1.5) abline(h=0,col=gray(0.3),lwd=1.5) par(mfrow=c(1,1)) n <- 9 pat <- rep("A",n) mat <- rep("a",n) pat mat mode(pat) from.mat <- rbinom(n,1,0.5) offspring <- pat from.mat==1 offspring[from.mat==1] mat[from.mat==1] offspring[from.mat==1] <- mat[from.mat==1] offspring 2:6 offspring[2:6] offspring[-(2:6)] meiosis <- function(GF,GM) { from.GM <- rbinom(length(GF),1,0.5) GS <- GF GS[from.GM==1] <- GM[from.GM==1] return(GS) } meiosis(pat,mat) n <- 10^5 model <- list(mu=5, alpha=1,delta=-1,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) 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 y <- m + model$sigma*rnorm(length(x)) return(list(pat=pat, mat=mat, pheno=y)) } F1 <- cross(IB1,IB2,model) BC1 <- cross(IB1,F1,model) BC2 <- cross(IB2,F1,model) F2 <- cross(F1,F1,model) 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) } op <- par(mfrow=c(2,3)) plot.cross("IB1",IB1) plot.cross("IB2",IB2) plot.cross("F1",F1) plot.cross("BC1",BC1) plot.cross("BC2",BC2) plot.cross("F2",F2) par(op)