load("nochr20.RData") library(dplyr) library(tidyr) library(ggplot2) library(ggpubr) 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) #observed value share_real <- vector(length=length(simdata[,1])) rgrm <- vector(length=length(simdata[,1])) f1 <- vector(length=length(simdata[,1])) f2 <- vector(length=length(simdata[,1])) m1 <- vector(length=length(simdata[,1])) m2 <- vector(length=length(simdata[,1])) dam <- vector(length=length(simdata[,1])) sire <- vector(length=length(simdata[,1])) rdmin <- vector(length=length(simdata[,1])) rdsum <- vector(length=length(simdata[,1])) rdmax <- vector(length=length(simdata[,1])) r0 <- 0 r1 <- 0 r2 <- 0 rfa <- 0 rfb <- 0 rfc <- 0 rfd <- 0 rfe <- 0 rff <- 0 rfg <- 0 rfh <- 0 rma <- 0 rmb <- 0 rmc <- 0 rmd <- 0 rme <- 0 rmf <- 0 rmg <- 0 rmh <- 0 rgrm_row <- 0 rgrm_zero <- 0 rhef <- 0 rhem <- 0 hof <- 0 hef <- 0 hom <- 0 hem <- 0 for (n in 1:length(simdata[,1])){ print(n) s <- 0 homo <- 0 dist1 <- 0 dist2 <- 0 dist3 <- 0 dist4 <- 0 f1[n] <- as.character(simdata[n,"f1"]) f2[n] <- as.character(simdata[n,"f2"]) m1[n] <- as.character(simdata[n,"m1"]) m2[n] <- as.character(simdata[n,"m2"]) dam[n] <- as.character(simdata[n,"dam"]) sire[n] <- as.character(simdata[n,"sire"]) rgrm_row <- which((grm$row==dam[n])&(grm$col==sire[n])) dist1 <- div$dist[which((div$h1==f1[n])&(div$h2==m1[n]))] dist2 <- div$dist[which((div$h1==f1[n])&(div$h2==m2[n]))] dist3 <- div$dist[which((div$h1==f2[n])&(div$h2==m1[n]))] dist4 <- div$dist[which((div$h1==f2[n])&(div$h2==m2[n]))] rdmin[n] <- min(dist1,dist2,dist3,dist4) rdsum[n] <- sum(dist1,dist2,dist3,dist4) rdmax[n] <- max(dist1,dist2,dist3,dist4) if (length(rgrm_row)==0){ rgrm_row <- which((grm$row==sire[n])&(grm$col==dam[n])) } if (length(rgrm_row)==0){ rgrm[n] <- 0 rgrm_zero <- rgrm_zero + 1 } else { rgrm[n] <- grm$number[rgrm_row] print (rgrm[n]) } if (f1[n]==f2[n]) { homo <- homo + 1 hof <- hof + 1 } else { homo <- homo +0 hef <- hef + 1} if (m1[n]==m2[n]) { homo <- homo +1 hom <- hom +1 } else { homo <- homo +0 hem <- hem +1} while (homo==0) { if (f1[n]==m1[n]){ s <- s+1 } if (f1[n]==m2[n]){ s <- s+1 } if (f2[n]==m1[n]){ s<- s+1 } if (f2[n]==m2[n]){ s<- s+1 } homo <- 3 } while (homo==1){ if (f1[n]==f2[n]){ if (f1[n]==m1[n]){ s <- s+1 } if (f1[n]==m2[n]){ s<- s+1 } } else { if (f1[n]==m1[n]){ s <- s+1 } if (f2[n]==m1[n]){ s<- s+1 } } homo <- 3 } while (homo==2){ if (f1[n]==m1[n]){ s <- 2 } homo <- 3 } share_real[n] <- s if (share_real[n]==0){ r0 <- r0 + 1 } if (share_real[n]==1){ r1 <- r1 + 1 } if (share_real[n]==2){ r2 <- r2 + 1 } if (f1[n]=="A") { rfa <- rfa + 1 } if (f2[n]=="A") { rfa <- rfa + 1 } if (f1[n]=="B") { rfb <- rfb + 1 } if (f2[n]=="B") { rfb <- rfb + 1 } if (f1[n]=="C") { rfc <- rfc + 1 } if (f2[n]=="C") { rfc <- rfc + 1 } if (f1[n]=="D") { rfd <- rfd + 1 } if (f2[n]=="D") { rfd <- rfd + 1 } if (f1[n]=="E") { rfe <- rfe + 1 } if (f2[n]=="E") { rfe <- rfe + 1 } if (f1[n]=="F") { rff <- rff + 1 } if (f2[n]=="F") { rff <- rff + 1 } if (f1[n]=="G") { rfg <- rfg + 1 } if (f2[n]=="G") { rfg <- rfg + 1 } if (f1[n]=="H") { rfh <- rfh + 1 } if (f2[n]=="H") { rfh <- rfh + 1 } if (m1[n]=="A") { rma <- rma + 1 } if (m2[n]=="A") { rma <- rma + 1 } if (m1[n]=="B") { rmb <- rmb + 1 } if (m2[n]=="B") { rmb <- rmb + 1 } if (m1[n]=="C") { rmc <- rmc + 1 } if (m2[n]=="C") { rmc <- rmc + 1 } if (m1[n]=="D") { rmd <- rmd + 1 } if (m2[n]=="D") { rmd <- rmd + 1 } if (m1[n]=="E") { rme <- rme + 1 } if (m2[n]=="E") { rme <- rme + 1 } if (m1[n]=="F") { rmf <- rmf + 1 } if (m2[n]=="F") { rmf <- rmf + 1 } if (m1[n]=="G") { rmg <- rmg + 1 } if (m2[n]=="G") { rmg <- rmg + 1 } if (m1[n]=="H") { rmh <- rmh + 1 } if (m2[n]=="H") { rmh <- rmh + 1 } } average_real <- sum(share_real)/(length(simdata[,1])) rgrm_real <- sum(rgrm)/(length(simdata[,1])) rgrm_ratio <- (rgrm_zero)/(length(simdata[,1])) simdata$grm <- rgrm rgrm_median <- median(rgrm) r0 <- r0/length(simdata[,1]) r1 <- r1/length(simdata[,1]) r2 <- r2/length(simdata[,1]) rfa <- rfa/(length(simdata[,1])*2) rfb <- rfb/(length(simdata[,1])*2) rfc <- rfc/(length(simdata[,1])*2) rfd <- rfd/(length(simdata[,1])*2) rfe <- rfe/(length(simdata[,1])*2) rff <- rff/(length(simdata[,1])*2) rfg <- rfg/(length(simdata[,1])*2) rfh <- rfh/(length(simdata[,1])*2) rma <- rma/(length(simdata[,1])*2) rmb <- rmb/(length(simdata[,1])*2) rmc <- rmc/(length(simdata[,1])*2) rmd <- rmd/(length(simdata[,1])*2) rme <- rme/(length(simdata[,1])*2) rmf <- rmf/(length(simdata[,1])*2) rmg <- rmg/(length(simdata[,1])*2) rmh <- rmh/(length(simdata[,1])*2) rhef <- hof/hef rhem <- hom/hem rmeandmin <- mean(rdmin) rmeandsum <- mean(rdsum) rmeandmax <- mean(rdmax) rmediandmin <- median(rdmin) rmediandsum <- median(rdsum) rmediandmax <- mean(rdmax) #choose the simulated result result <- read.csv("adjust_simu.csv") result <- read.csv("random_simu.csv") result <- read.csv("restrict_simu.csv") 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 consort1 <- read.csv("consort1.csv",header=T) sim1 <- consort1 share_real1 <- vector(length=length(sim1[,1])) rgrm1 <- vector(length=length(sim1[,1])) f1 <- vector(length=length(sim1[,1])) f2 <- vector(length=length(sim1[,1])) m1 <- vector(length=length(sim1[,1])) m2 <- vector(length=length(sim1[,1])) dam <- vector(length=length(sim1[,1])) sire <- vector(length=length(sim1[,1])) rdmin1 <- vector(length=length(simdata[,1])) rdsum1<- vector(length=length(simdata[,1])) rdmax1<- vector(length=length(simdata[,1])) r0a <- 0 r1a <- 0 r2a <- 0 rfa1 <- 0 rfb1 <- 0 rfc1 <- 0 rfd1 <- 0 rfe1 <- 0 rff1 <- 0 rfg1 <- 0 rfh1 <- 0 rma1 <- 0 rmb1 <- 0 rmc1 <- 0 rmd1 <- 0 rme1 <- 0 rmf1 <- 0 rmg1 <- 0 rmh1 <- 0 rgrm_row1 <- 0 rgrm1 <- 0 rgrm_zero1 <- 0 rhef1 <- 0 rhem1 <- 0 hof1 <- 0 hef1 <- 0 hom1 <- 0 hem1 <- 0 for (n in 1:length(sim1[,1])){ print(n) s <- 0 homo <- 0 dist1 <- 0 dist2 <- 0 dist3 <- 0 dist4 <- 0 f1[n] <- as.character(sim1[n,"f1"]) f2[n] <- as.character(sim1[n,"f2"]) m1[n] <- as.character(sim1[n,"m1"]) m2[n] <- as.character(sim1[n,"m2"]) dam[n] <- as.character(sim1[n,"mother"]) sire[n] <- as.character(sim1[n,"father"]) dist1 <- div$dist[which((div$h1==f1[n])&(div$h2==m1[n]))] dist2 <- div$dist[which((div$h1==f1[n])&(div$h2==m2[n]))] dist3 <- div$dist[which((div$h1==f2[n])&(div$h2==m1[n]))] dist4 <- div$dist[which((div$h1==f2[n])&(div$h2==m2[n]))] rdmin1[n] <- min(dist1,dist2,dist3,dist4) rdsum1[n] <- sum(dist1,dist2,dist3,dist4) rdmax1[n] <- max(dist1,dist2,dist3,dist4) if (f1[n]==f2[n]) { homo <- homo + 1 hof1 <- hof1 + 1 } else { homo <- homo +0 hef1 <- hef1 + 1} if (m1[n]==m2[n]) { homo <- homo +1 hom1 <- hom1 +1 } else { homo <- homo +0 hem1 <- hem1 +1} while (homo==0) { if (f1[n]==m1[n]){ s <- s+1 } if (f1[n]==m2[n]){ s <- s+1 } if (f2[n]==m1[n]){ s<- s+1 } if (f2[n]==m2[n]){ s<- s+1 } homo <- 3 } while (homo==1){ if (f1[n]==f2[n]){ if (f1[n]==m1[n]){ s <- s+1 } if (f1[n]==m2[n]){ s<- s+1 } } else { if (f1[n]==m1[n]){ s <- s+1 } if (f2[n]==m1[n]){ s<- s+1 } } homo <- 3 } while (homo==2){ if (f1[n]==m1[n]){ s <- 2 } homo <- 3 } share_real1[n] <- s if (share_real1[n]==0){ r0a <- r0a + 1 } if (share_real1[n]==1){ r1a <- r1a + 1 } if (share_real1[n]==2){ r2a <- r2a + 1 } if (f1[n]=="A") { rfa1 <- rfa1 + 1 } if (f2[n]=="A") { rfa1 <- rfa1 + 1 } if (f1[n]=="B") { rfb1 <- rfb1 + 1 } if (f2[n]=="B") { rfb1 <- rfb1 + 1 } if (f1[n]=="C") { rfc1 <- rfc1 + 1 } if (f2[n]=="C") { rfc1 <- rfc1 + 1 } if (f1[n]=="D") { rfd1 <- rfd1 + 1 } if (f2[n]=="D") { rfd1 <- rfd1 + 1 } if (f1[n]=="E") { rfe1 <- rfe1 + 1 } if (f2[n]=="E") { rfe1 <- rfe1 + 1 } if (f1[n]=="F") { rff1 <- rff1 + 1 } if (f2[n]=="F") { rff1 <- rff1 + 1 } if (f1[n]=="G") { rfg1 <- rfg1 + 1 } if (f2[n]=="G") { rfg1 <- rfg1 + 1 } if (f1[n]=="H") { rfh1 <- rfh1 + 1 } if (f2[n]=="H") { rfh1 <- rfh1 + 1 } if (m1[n]=="A") { rma1 <- rma1 + 1 } if (m2[n]=="A") { rma1 <- rma1 + 1 } if (m1[n]=="B") { rmb1 <- rmb1 + 1 } if (m2[n]=="B") { rmb1 <- rmb1 + 1 } if (m1[n]=="C") { rmc1 <- rmc1 + 1 } if (m2[n]=="C") { rmc1 <- rmc1 + 1 } if (m1[n]=="D") { rmd1 <- rmd1 + 1 } if (m2[n]=="D") { rmd1<- rmd1 + 1 } if (m1[n]=="E") { rme1 <- rme1 + 1 } if (m2[n]=="E") { rme1 <- rme1 + 1 } if (m1[n]=="F") { rmf1 <- rmf1 + 1 } if (m2[n]=="F") { rmf1 <- rmf1 + 1 } if (m1[n]=="G") { rmg1 <- rmg1 + 1 } if (m2[n]=="G") { rmg1 <- rmg1 + 1 } if (m1[n]=="H") { rmh1 <- rmh1 + 1 } if (m2[n]=="H") { rmh1 <- rmh1 + 1 } } average_real1 <- sum(share_real1)/(length(sim1[,1])) rgrm_real1 <- sum(sim1$grm)/(length(sim1[,1])) rgrm_median1 <- median(sim1$grm) r0a <- r0a/length(sim1[,1]) r1a <- r1a/length(sim1[,1]) r2a <- r2a/length(sim1[,1]) rfa1 <- rfa1/(length(sim1[,1])*2) rfb1 <- rfb1/(length(sim1[,1])*2) rfc1 <- rfc1/(length(sim1[,1])*2) rfd1 <- rfd1/(length(sim1[,1])*2) rfe1 <- rfe1/(length(sim1[,1])*2) rff1 <- rff1/(length(sim1[,1])*2) rfg1 <- rfg1/(length(sim1[,1])*2) rfh1 <- rfh1/(length(sim1[,1])*2) rma1 <- rma1/(length(sim1[,1])*2) rmb1 <- rmb1/(length(sim1[,1])*2) rmc1 <- rmc1/(length(sim1[,1])*2) rmd1 <- rmd1/(length(sim1[,1])*2) rme1 <- rme1/(length(sim1[,1])*2) rmf1 <- rmf1/(length(sim1[,1])*2) rmg1 <- rmg1/(length(sim1[,1])*2) rmh1 <- rmh1/(length(sim1[,1])*2) rhef1 <- hof1/hef1 rhem1 <- hom1/hem1 rmeandmin1 <- mean(rdmin1) rmeandsum1 <- mean(rdsum1) rmeandmax1 <- mean(rdmax1) rmediandmin1 <- median(rdmin1) rmediandsum1 <- median(rdsum1) rmediandmax1 <- median(rdmax1) data <- data.frame(average) data1 <- data.frame(a) data2 <- data.frame(b) data3 <- data.frame(c) data4 <- data.frame(fga,fgb,fgc,fgd,fge,fgf,fgg,fgh) data5 <- data.frame(mga,mgb,mgc,mgd,mge,mgf,mgg,mgh) data6 <- data.frame(grm_average) data7 <- data.frame(grm_median) data8 <- data.frame(hehof,hehom) data9 <- data.frame(meandmax,meandsum,mediandmax,mediandsum) p1 <- ggplot(aes(x=average), data=data) + geom_vline(xintercept=average_real,lty=2,color="blue")+ geom_vline(xintercept=average_real1,lty=2,color="red") + geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Average number of shared haplotype per mating pair")+geom_vline(xintercept=quantile(average,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(average,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())+ggtitle("c") p2 <- ggplot(aes(x=a), data=data1) + geom_vline(xintercept=r0,lty=2,color="blue") + geom_vline(xintercept=r0a,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Proportion of mating pairs sharing 0 haplotype")+geom_vline(xintercept=quantile(a,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(a,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())+ggtitle("d") p3 <- ggplot(aes(x=b), data=data2) + geom_vline(xintercept=r1,lty=2,color="blue") + geom_vline(xintercept=r1a,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Proportion of mating pairs sharing 1 haplotype")+geom_vline(xintercept=quantile(b,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(b,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())+ggtitle("e") p4 <- ggplot(aes(x=c), data=data3) + geom_vline(xintercept=r2,lty=2,color="blue") + geom_vline(xintercept=r2a,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Proportion of mating pairs sharing 2 haplotypes")+geom_vline(xintercept=quantile(c,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(c,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())+ggtitle("f") p5 <- ggplot(aes(x=grm_average), data=data6) + geom_vline(xintercept=rgrm_real1,lty=2,color="red")+ geom_vline(xintercept=rgrm_real,lty=2,color="blue")+geom_histogram(binwidth = 0.0001,alpha=0.8) + scale_x_continuous(breaks = seq(-1, 1, by = 0.002))+xlab("Average of genomic relatedness across all mating pairs")+geom_vline(xintercept=quantile(grm_average,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(grm_average,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())+ggtitle("a") p6 <- ggplot(aes(x=grm_median), data=data7) + geom_vline(xintercept=rgrm_median1,lty=2,color="red") + geom_vline(xintercept=rgrm_median,lty=2,color="blue")+geom_histogram(binwidth = 0.0001,alpha=0.8) + scale_x_continuous(breaks = seq(-0.02, 0, by = 0.002))+xlab("Median of genomic relatedness across all mating pairs")+geom_vline(xintercept=quantile(grm_median,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(grm_median,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())+ggtitle("b") p7 <- ggplot(aes(x=hehof), data=data8) + geom_vline(xintercept=rhef,lty=2,color="blue")+ geom_vline(xintercept=rhef1,lty=2,color="red")+ geom_histogram(binwidth = 0.002,alpha=0.8) + scale_x_continuous(breaks = seq(0, 0.19, by = 0.01))+xlab("Ratio of homozygote/heterozygote in mothers")+geom_vline(xintercept=quantile(hehof,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(hehof,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())+ggtitle("a") p8 <- ggplot(aes(x=hehom), data=data8) + geom_vline(xintercept=rhem,lty=2,color="blue")+ geom_vline(xintercept=rhem1,lty=2,color="red")+ geom_histogram(binwidth = 0.002,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Ratio of homozygot/heterozygote in fathers")+geom_vline(xintercept=quantile(hehom,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(hehom,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())+ggtitle("b") p9 <- ggplot(aes(x=meandsum), data=data9) + geom_vline(xintercept=rmeandsum,lty=2,color="blue")+ geom_vline(xintercept=rmeandsum1,lty=2,color="red")+ geom_histogram(binwidth = 0.0002,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.005))+xlab("Mean of AAdist across all mating pairs")+geom_vline(xintercept=quantile(meandsum,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(meandsum,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())+ggtitle("g") p10 <- ggplot(aes(x=meandmax), data=data9) + geom_vline(xintercept=rmeandmax,lty=2,color="blue")+ geom_vline(xintercept=rmeandmax1,lty=2,color="red")+ geom_histogram(binwidth = 0.0002,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.002))+xlab("Mean of distmax across all mating pairs")+geom_vline(xintercept=quantile(meandmax,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(meandmax,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())+ggtitle("h") fa <- ggplot(aes(x=fga), data=data4) + geom_vline(xintercept=rfa,lty=2,color="blue")+ geom_vline(xintercept=rfa1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype A in mothers")+ggtitle("A")+geom_vline(xintercept=quantile(fga,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fga,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) fb <- ggplot(aes(x=fgb), data=data4) + geom_vline(xintercept=rfb,lty=2,color="blue")+ geom_vline(xintercept=rfb1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype B in mothers")+ggtitle("B")+geom_vline(xintercept=quantile(fgb,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fgb,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) fc <- ggplot(aes(x=fgc), data=data4) + geom_vline(xintercept=rfc,lty=2,color="blue")+ geom_vline(xintercept=rfc1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype C in mothers ")+ggtitle("C")+geom_vline(xintercept=quantile(fgc,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fgc,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) fd <- ggplot(aes(x=fgd), data=data4) + geom_vline(xintercept=rfd,lty=2,color="blue")+ geom_vline(xintercept=rfd1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype D in mothers")+ggtitle("D")+geom_vline(xintercept=quantile(fgd,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fgd,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) fe <- ggplot(aes(x=fge), data=data4) + geom_vline(xintercept=rfe,lty=2,color="blue")+ geom_vline(xintercept=rfe1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype E in mothers")+ggtitle("E")+geom_vline(xintercept=quantile(fge,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fge,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) ff <- ggplot(aes(x=fgf), data=data4) + geom_vline(xintercept=rff,lty=2,color="blue")+ geom_vline(xintercept=rff1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype F in mothers")+ggtitle("F")+geom_vline(xintercept=quantile(fgf,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fgf,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) fg <- ggplot(aes(x=fgg), data=data4) + geom_vline(xintercept=rfg,lty=2,color="blue")+ geom_vline(xintercept=rfg1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype G in mothers")+ggtitle("G")+geom_vline(xintercept=quantile(fgg,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fgg,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) fh <- ggplot(aes(x=fgh), data=data4) + geom_vline(xintercept=rfh,lty=2,color="blue")+ geom_vline(xintercept=rfh1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype H in mothers")+ggtitle("H")+geom_vline(xintercept=quantile(fgh,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fgh,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) ma <- ggplot(aes(x=mga), data=data5) + geom_vline(xintercept=rma,lty=2,color="blue")+ geom_vline(xintercept=rma1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype A in fathers")+ggtitle("A")+geom_vline(xintercept=quantile(mga,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mga,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) mb <- ggplot(aes(x=mgb), data=data5) + geom_vline(xintercept=rmb,lty=2,color="blue")+ geom_vline(xintercept=rmb1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype B in fathers")+ggtitle("B")+geom_vline(xintercept=quantile(mgb,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mgb,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) mc <- ggplot(aes(x=mgc), data=data5) + geom_vline(xintercept=rmc,lty=2,color="blue")+ geom_vline(xintercept=rmc1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype C in fathers")+ggtitle("C")+geom_vline(xintercept=quantile(mgc,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mgc,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) md <- ggplot(aes(x=mgd), data=data5) + geom_vline(xintercept=rmd,lty=2,color="blue")+ geom_vline(xintercept=rmd1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype D in fathers")+ggtitle("D")+geom_vline(xintercept=quantile(mgd,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mgd,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) me <- ggplot(aes(x=mge), data=data5) + geom_vline(xintercept=rme,lty=2,color="blue")+ geom_vline(xintercept=rme1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype E in fathers")+ggtitle("E")+geom_vline(xintercept=quantile(mge,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mge,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) mf <- ggplot(aes(x=mgf), data=data5) + geom_vline(xintercept=rmf,lty=2,color="blue")+ geom_vline(xintercept=rmf1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype F in fathers")+ggtitle("F")+geom_vline(xintercept=quantile(mgf,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mgf,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) mg <- ggplot(aes(x=mgg), data=data5) + geom_vline(xintercept=rmg,lty=2,color="blue")+ geom_vline(xintercept=rmg1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype G in fathers")+ggtitle("G")+geom_vline(xintercept=quantile(mgg,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mgg,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) mh <- ggplot(aes(x=mgh), data=data5) + geom_vline(xintercept=rmh,lty=2,color="blue")+ geom_vline(xintercept=rmh1,lty=2,color="red")+ geom_histogram(binwidth = 0.001,alpha=0.8) + scale_x_continuous(breaks = seq(0, 1, by = 0.02))+xlab("Frequency of haplotype H in fathers")+ggtitle("H")+geom_vline(xintercept=quantile(mgh,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mgh,0.9984375),lty=3,size=1)+theme(axis.title.y = element_blank()) p11 <- ggarrange(fa,fb,fc,fd,fe,ff,fg,fh,nrow=4,ncol=2) p12 <- ggarrange(ma,mb,mc,md,me,mf,mg,mh,nrow=4,ncol=2) annotate_figure(p11, left="Number of simulations") annotate_figure(p12, left="Number of simulations") p13 <- ggarrange(p7,p8,p1,p2,p3,p4,p9,p10,nrow=4,ncol=2) annotate_figure(p13, left="Number of simulations") p14 <- ggarrange(p5,p6,nrow=1,ncol=2) annotate_figure(p14, left="Number of simulations")