# Inbreeding at Several Loci 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,rec.frac) { pat <- meiosis.rec(fa$pat,fa$mat,rec.frac) mat <- meiosis.rec(mo$pat,mo$mat,rec.frac) return(list(pat=pat,mat=mat)) } 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)) IB2 <- list(pat=cbind(A,B),mat=cbind(A,B)) F1 <- cross.rec(IB1,IB2,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,theta) theta.BC <- c(theta.BC,rec.count(BC2)) F2 <- cross.rec(F1,F1,theta) theta.F2 <- c(theta.F2,rec.count(F2)) Fn.fa <- cross.rec(F1,F1,theta) Fn.mo <- cross.rec(F1,F1,theta) for (g in 2:20) { New.fa <- cross.rec(Fn.fa,Fn.mo,theta) New.mo <- cross.rec(Fn.fa,Fn.mo,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)) theta.RI <- c(theta.RI, rec.count(RI)) } plot(rec.frac,theta.RI,type="l",xlab="theta", ylab="rec. fraction") lines(rec.frac,theta.F2,col=2) lines(rec.frac,theta.BC,col=3) legend(0,0.5,legend=c("RI","F2","BC"),lty=rep(1,3),col=1:3) lines(rec.frac,4*rec.frac/(1+6*rec.frac),lty=2)