load("nochr20.rData") library(dplyr) library(tidyr) simdata <- read.csv("offspring_NF.csv",header=T) male <- read.csv("adult_male_NF.csv",header=T) female <- read.csv("adult_female_NF.csv",header=T) div <- read.csv("MHC.csv") grm <- NoChr20GRM_030420 colnames(grm) <- c("row","col","number") grm3 <- subset(grm,grm$row%in%male$sire & grm$col%in%female$dam) grm4 <- subset(grm,grm$row%in%female$dam & grm$col%in%male$sire) grm <- rbind(grm3,grm4) sexual_selection <- function(nsim){ average <- vector(length=nsim) a <- vector(length=nsim) b <- vector(length=nsim) c <- vector(length=nsim) fga <- vector(length=nsim) fgb <- vector(length=nsim) fgc <- vector(length=nsim) fgd <- vector(length=nsim) fge <- vector(length=nsim) fgf <- vector(length=nsim) fgg <- vector(length=nsim) fgh <- vector(length=nsim) mga <- vector(length=nsim) mgb <- vector(length=nsim) mgc <- vector(length=nsim) mgd <- vector(length=nsim) mge <- vector(length=nsim) mgf <- vector(length=nsim) mgg <- vector(length=nsim) mgh <- vector(length=nsim) grm_total <- vector(length=nsim) grm_ratio <- vector(length=nsim) grm_median <- vector(length=nsim) grm_average <- vector(length=nsim) sum_mabs <- NULL sum_fabs <- NULL sum_mlbs <- NULL sum_flbs <- NULL hehof <- vector(length=nsim) hehom <- vector(length=nsim) meandmin <- vector(length=nsim) meandsum <- vector(length=nsim) meandmax <- vector(length=nsim) mediandmin <- vector(length=nsim) mediandsum <- vector(length=nsim) mediandmax <- vector(length=nsim) for (j in 1:nsim){ simu1 <- NULL mabs <- NULL fabs <- NULL mlbs <- NULL flbs <- NULL print(j) for (m in 1988:2011){ male1 <- male[male$year==m,] female1 <- female[female$year==m,] o <- m + 1 simu <- simdata[simdata$birthyear==o,] #sample from mating pool x <- sample(1:length(male1[,1]),length(simu[,1]),replace=T) y <- sample(1:length(female1[,1]),length(simu[,1]),replace=F) fatherid <- vector(length=length(simu[,1])) motherid <- vector(length=length(simu[,1])) h1 <- vector(length=length(simu[,1])) h2 <- vector(length=length(simu[,1])) h3 <- vector(length=length(simu[,1])) h4 <- vector(length=length(simu[,1])) animal <- vector(length=length(simu[,1])) grm_simu <- vector(length=length(simu[,1])) simu2 <- NULL mbs <- NULL fbs <- NULL print(m) for (i in 1:length(simu[,1])){ fatherid[i] <- male1$sire[x[i]] motherid[i] <- female1$dam[y[i]] h1[i] <- as.character(male1$h1[x[i]]) h2[i] <- as.character(male1$h2[x[i]]) h3[i] <- as.character(female1$h1[y[i]]) h4[i] <- as.character(female1$h2[y[i]]) animal[i] <- simu$animal[i] grm_row <- which((grm$row==fatherid[i])&(grm$col==motherid[i])) if (length(grm_row)==0) {grm_row <- which((grm$col==fatherid[i])&(grm$row==motherid[i]))} if (length(grm_row)==0){ grm_simu[i] <- 0 } else { grm_simu[i] <- grm$number[grm_row] } } male_sample <- data.frame(animal,fatherid,h1,h2,grm_simu) female_sample <- data.frame(animal,motherid,h3,h4) simu2 <- left_join(female_sample,male_sample,by="animal") simu2$yr <- rep(m,length(simu[,1])) simu1 <- rbind(simu1,simu2) mbs<- as.data.frame(table(simu2$fatherid)) fbs<- as.data.frame(table(simu2$motherid)) mbs$yr <- rep(m,length(mbs[,1])) fbs$yr <- rep(m,length(fbs[,1])) mabs <- rbind(mabs,mbs) fabs <- rbind(fabs,fbs) } simu1$cycle <- rep(j,length(simdata[,1])) mlbs<- as.data.frame(table(simu1$fatherid)) flbs<- as.data.frame(table(simu1$motherid)) flbs$cycle <- rep(j,length(flbs[,1])) mlbs$cycle <- rep(j,length(mlbs[,1])) fabs$cycle <- rep(j,length(fabs[,1])) mabs$cycle <- rep(j,length(mabs[,1])) sum_fabs <- rbind(sum_fabs,fabs) sum_mabs <- rbind(sum_mabs,mabs) sum_flbs <- rbind(sum_flbs,flbs) sum_mlbs <- rbind(sum_mlbs,mlbs) share <- vector(length=length(simu1[,1])) f1 <- vector(length=length(simu1[,1])) f2 <- vector(length=length(simu1[,1])) m1 <- vector(length=length(simu1[,1])) m2 <- vector(length=length(simu1[,1])) grm_sum <- vector(length=length(simu1[,1])) dmin <- vector(length=length(simu1[,1])) dsum <- vector(length=length(simu1[,1])) dmax <- vector(length=length(simu1[,1])) s0 <- 0 s1 <- 0 s2 <- 0 fa <- 0 fb <- 0 fc <- 0 fd <- 0 fe <- 0 ff <- 0 fg <- 0 fh <- 0 ma <- 0 mb <- 0 mc <- 0 md <- 0 me <- 0 mf <- 0 mg <- 0 mh <- 0 grm_p <- 0 hetef <- 0 hetem <- 0 homof <- 0 homom <- 0 for (k in 1:length(simu1[,1])) { s <- 0 homo <- 0 dist1 <- 0 dist2 <- 0 dist3 <- 0 dist4 <- 0 f1[k] <- as.character(simu1[k,"h3"]) f2[k] <- as.character(simu1[k,"h4"]) m1[k] <- as.character(simu1[k,"h1"]) m2[k] <- as.character(simu1[k,"h2"]) grm_sum[k] <- as.numeric(simu1[k,"grm_simu"]) dist1 <- div$dist[which((div$h1==f1[k])&(div$h2==m1[k]))] dist2 <- div$dist[which((div$h1==f1[k])&(div$h2==m2[k]))] dist3 <- div$dist[which((div$h1==f2[k])&(div$h2==m1[k]))] dist4 <- div$dist[which((div$h1==f2[k])&(div$h2==m2[k]))] dmin[k] <- min(dist1,dist2,dist3,dist4) dsum[k] <- sum(dist1,dist2,dist3,dist4) dmax[k] <- max(dist1,dist2,dist3,dist4) if (grm_sum[k]==0) { grm_p <- grm_p +1 } if (f1[k]==f2[k]) { homo <- homo + 1 homof <- homof +1 } else { homo <- homo +0 hetef <- hetef +1} if (m1[k]==m2[k]) { homo <- homo +1 homom <- homom +1 } else { homo <- homo +0 hetem <- hetem +1} while (homo==0) { if (f1[k]==m1[k]){ s <- s+1 } if (f1[k]==m2[k]){ s <- s+1 } if (f2[k]==m1[k]){ s<- s+1 } if (f2[k]==m2[k]){ s<- s+1 } homo <- 3 } while (homo==1){ if (f1[k]==f2[k]){ if (f1[k]==m1[k]){ s <- s+1 } if (f1[k]==m2[k]){ s<- s+1 } } else { if (f1[k]==m1[k]){ s <- s+1 } if (f2[k]==m1[k]){ s<- s+1 } } homo <- 3 } while (homo==2){ if (f1[k]==m1[k]){ s <- 2 } homo <- 3 } share[k] <- s if (share[k]==0){ s0 <- s0 + 1 } if (share[k]==1){ s1 <- s1 + 1 } if (share[k]==2){ s2 <- s2 + 1 } if (f1[k]=="A") { fa <- fa + 1 } if (f2[k]=="A") { fa <- fa + 1 } if (f1[k]=="B") { fb <- fb + 1 } if (f2[k]=="B") { fb <- fb + 1 } if (f1[k]=="C") { fc <- fc + 1 } if (f2[k]=="C") { fc <- fc + 1 } if (f1[k]=="D") { fd <- fd + 1 } if (f2[k]=="D") { fd <- fd + 1 } if (f1[k]=="E") { fe <- fe + 1 } if (f2[k]=="E") { fe <- fe + 1 } if (f1[k]=="F") { ff <- ff + 1 } if (f2[k]=="F") { ff <- ff + 1 } if (f1[k]=="G") { fg <- fg + 1 } if (f2[k]=="G") { fg <- fg + 1 } if (f1[k]=="H") { fh <- fh + 1 } if (f2[k]=="H") { fh <- fh + 1 } if (m1[k]=="A") { ma <- ma + 1 } if (m2[k]=="A") { ma <- ma + 1 } if (m1[k]=="B") { mb <- mb + 1 } if (m2[k]=="B") { mb <- mb + 1 } if (m1[k]=="C") { mc <- mc + 1 } if (m2[k]=="C") { mc <- mc + 1 } if (m1[k]=="D") { md <- md + 1 } if (m2[k]=="D") { md <- md + 1 } if (m1[k]=="E") { me <- me + 1 } if (m2[k]=="E") { me <- me + 1 } if (m1[k]=="F") { mf <- mf + 1 } if (m2[k]=="F") { mf <- mf + 1 } if (m1[k]=="G") { mg <- mg + 1 } if (m2[k]=="G") { mg <- mg + 1 } if (m1[k]=="H") { mh <- mh + 1 } if (m2[k]=="H") { mh <- mh + 1 } } average[j] <- sum(share)/(length(simu1[,1])) a[j] <- s0/(length(simu1[,1])) b[j] <- s1/(length(simu1[,1])) c[j] <- s2/(length(simu1[,1])) fga[j] <- fa/((length(simu1[,1]))*2) fgb[j] <- fb/((length(simu1[,1]))*2) fgc[j] <- fc/((length(simu1[,1]))*2) fgd[j] <- fd/((length(simu1[,1]))*2) fge[j] <- fe/((length(simu1[,1]))*2) fgf[j] <- ff/((length(simu1[,1]))*2) fgg[j] <- fg/((length(simu1[,1]))*2) fgh[j] <- fh/((length(simu1[,1]))*2) mga[j] <- ma/((length(simu1[,1]))*2) mgb[j] <- mb/((length(simu1[,1]))*2) mgc[j] <- mc/((length(simu1[,1]))*2) mgd[j] <- md/((length(simu1[,1]))*2) mge[j] <- me/((length(simu1[,1]))*2) mgf[j] <- mf/((length(simu1[,1]))*2) mgg[j] <- mg/((length(simu1[,1]))*2) mgh[j] <- mh/((length(simu1[,1]))*2) grm_total[j] <- sum(grm_sum) grm_ratio[j] <- (grm_p)/(length(simu1[,1])) grm_median[j] <- median(grm_sum) grm_average[j] <- grm_total[j]/(length(simu1[,1])) hehof[j] <- homof/hetef hehom[j] <- homom/hetem meandmin[j] <- mean(dmin) meandsum[j] <- mean(dsum) meandmax[j] <- mean(dmax) mediandmin[j] <- median(dmin) mediandsum[j] <- median(dsum) mediandmax[j] <- median(dmax) } mylist <- list("average"=average,"a"=a,"b"=b,"c"=c,"fga"=fga,"fgb"=fgb,"fgc"=fgc,"fgd"=fgd,"fge"=fge,"fgf"=fgf,"fgg"=fgg,"fgh"=fgh,"mga"=mga,"mgb"=mgb,"mgc"=mgc,"mgd"=mgd,"mge"=mge,"mgf"=mgf,"mgg"=mgg,"mgh"=mgh,"grm_median"=grm_median,"grm_average"=grm_average,"hehof"=hehof,"hehom"=hehom,"mediandmin"=mediandmin,"mediandsum"=mediandsum,"mediandmax"=mediandmax,"meandmin"=meandmin,"meandsum"=meandsum,"meandmax"=meandmax,"sum_fabs"=sum_fabs,"sum_mabs"=sum_mabs,"sum_flbs"=sum_flbs,"sum_mlbs"=sum_mlbs) } #run the simulation for 10000 iterations result <- sexual_selection(10000) average <- result$average a <- result$a b <- result$b c <- result$c fga <- result$fga fgb <- result$fgb fgc <- result$fgc fgd <- result$fgd fge <- result$fge fgf <- result$fgf fgg <- result$fgg fgh <- result$fgh mga <- result$mga mgb <- result$mgb mgc <- result$mgc mgd <- result$mgd mge <- result$mge mgf <- result$mgf mgg <- result$mgg mgh <- result$mgh grm_average <- result$grm_average grm_median <- result$grm_median hehof <- result$hehof hehom <- result$hehom meandmin <- result$meandmin mediandmin <- result$mediandmin meandsum <- result$meandsum mediandsum <- result$mediandsum meandmax <- result$meandmax mediandmax <- result$mediandmax sum_fabs <- result$sum_fabs sum_mabs <- result$sum_mabs sum_flbs <- result$sum_flbs sum_mlbs <- result$sum_mlbs result_summary <- data.frame(average,a,b,c,fga,fgb,fgc,fgd,fge,fgf,fgg,fgh,mga,mgb,mgc,mgd,mge,mgf,mgg,mgh,grm_median,grm_average,hehof,hehom,mediandmin,mediandsum,mediandmax,meandmin,meandsum,meandmax) sum_fabs1 <- filter(sum_fabs,sum_fabs$cycle < 6) sum_mabs1 <- filter(sum_mabs,sum_mabs$cycle < 6) sum_flbs1 <- filter(sum_flbs,sum_flbs$cycle < 6) sum_mlbs1 <- filter(sum_mlbs,sum_mlbs$cycle < 6) write.csv(result_summary,"random_simu.csv") write.csv(sum_fabs1,"fabs_random.csv") write.csv(sum_mabs1,"mabs_random.csv") write.csv(sum_flbs1,"flbs_random.csv") write.csv(sum_mlbs1,"mlbs_random.csv")