library(ggplot2)
library(ggpubr)
offspring <- read.csv("MHC_trio.csv",header=T)
div <- read.csv("MHC_divergence.csv")

#identify group based parental MHC diplotype combination 
f1 <- as.character(offspring$f1)
f2 <- as.character(offspring$f2)
m1 <- as.character(offspring$m1)
m2 <- as.character(offspring$m2)
type <- vector(length=length(offspring[,1]))
for (i in 1:length(offspring[,1])){
  print(i)
  
  if (f1[i]==f2[i] && m1[i]==m2[i]){
    if (f1[i]==m1[i]){
      type[i] <- "1" } else {
        type[i] <- "2"
      }
  }
  
  if (f1[i]==f2[i] && m1[i]!=m2[i]){
    if (f1[i]==m1[i] || f1[i]==m2[i]) {
      type[i]<- "5"
    } else{
      type[i]<- "3"
    }
  }
  
  if (f1[i]!=f2[i] && m1[i]==m2[i]){
    if (f1[i]==m1[i] || f2[i]==m1[i]) {
      type[i]<- "5"
    } else{
      type[i]<- "3"
    }
  }
  if (f1[i]!=f2[i] && m1[i]!=m2[i]){
    m <-0
    if (f1[i]==m1[i]) {m <- m +1}
    if (f1[i]==m2[i]) {m <- m +1}
    if (f2[i]==m1[i]) {m <- m +1}
    if (f2[i]==m2[i]) {m <- m +1}
    if (m == 0){type [i] <- "4"}
    if (m == 1){type [i] <- "7"} 
    if (m == 2){type [i] <- "6"}
  }
}

# remove group 1 and 2 for further analysis   
offspring <- offspring[offspring$type!="1"& offspring$type!="2",]

#Simulation
simulation_post <- function(nsim){
  # Analysis 1 : ratio of homozygote/heterozygote for each group 
  h5 <- vector(length= nsim)
  h6 <- vector(length= nsim)
  h7 <- vector(length= nsim)
  
  # Analysis 1 : number of heterozygote (he) and homozygote (ho) for each group 
  he31 <- vector(length= nsim)
  ho31 <- vector(length= nsim)
  he41 <- vector(length= nsim)
  ho41 <- vector(length= nsim)
  he51 <- vector(length= nsim)
  ho51 <- vector(length= nsim)
  he61 <- vector(length= nsim)
  ho61 <- vector(length= nsim)
  he71 <- vector(length= nsim)
  ho71 <- vector(length= nsim)
  
  # Analysis 2 : number of individuals with diplotype identical to father or mother in group 5 and 7 
  f51 <- vector(length= nsim)
  m51 <- vector(length= nsim)
  f71 <- vector(length= nsim)
  m71 <- vector(length= nsim)
  fm71 <- vector(length= nsim)
  
  # Analysis 4 (haplotype A-H)
  fga1 <- vector(length=nsim)
  fgb1 <- vector(length=nsim)
  fgc1 <- vector(length=nsim)
  fgd1 <- vector(length=nsim)
  fge1 <- vector(length=nsim)
  fgf1 <- vector(length=nsim)
  fgg1 <- vector(length=nsim)
  fgh1 <- vector(length=nsim)
  
  # Analysis 5 (number of each haplotype inherited from father (mht) or mother (fht))
  fhta1 <- vector(length=nsim)
  fhtb1 <- vector(length=nsim)
  fhtc1 <- vector(length=nsim)
  fhtd1 <- vector(length=nsim)
  fhte1 <- vector(length=nsim)
  fhtf1 <- vector(length=nsim)
  fhtg1 <- vector(length=nsim)
  fhth1 <- vector(length=nsim)
  mhta1 <- vector(length=nsim)
  mhtb1 <- vector(length=nsim)
  mhtc1 <- vector(length=nsim)
  mhtd1 <- vector(length=nsim)
  mhte1 <- vector(length=nsim)
  mhtf1 <- vector(length=nsim)
  mhtg1 <- vector(length=nsim)
  mhth1 <- vector(length=nsim)
  
  # Analysis 6 (number of each haplotype received in males (omh) or females (ofh))
  ofha1 <- vector(length=nsim)
  ofhb1 <- vector(length=nsim)
  ofhc1 <- vector(length=nsim)
  ofhd1 <- vector(length=nsim)
  ofhe1 <- vector(length=nsim)
  ofhf1 <- vector(length=nsim)
  ofhg1 <- vector(length=nsim)
  ofhh1 <- vector(length=nsim)
  omha1 <- vector(length=nsim)
  omhb1 <- vector(length=nsim)
  omhc1 <- vector(length=nsim)
  omhd1 <- vector(length=nsim)
  omhe1 <- vector(length=nsim)
  omhf1 <- vector(length=nsim)
  omhg1 <- vector(length=nsim)
  omhh1 <- vector(length=nsim)
  
  # Analysis 3 mean (dist) and median (distm) MHC divergence
  dist <- vector(length=nsim)
  distm <- vector(length=nsim)
  
  for (j in 1: nsim){
    print(j)
    
    # simulated MHC haplotype from mother (o) and father (p)  
    o <- vector(length = length(offspring[,1]))
    p <- vector(length = length(offspring[,1]))
    
    # real MHC haplotype of father (m1 and m2) and mother (f1 and f2) 
    fz1 <- vector(length = length(offspring[,1]))
    mz1 <- vector(length = length(offspring[,1]))
    f1 <- vector(length = length(offspring[,1]))
    m1 <- vector(length = length(offspring[,1]))
    f2 <- vector(length = length(offspring[,1]))
    m2 <- vector(length = length(offspring[,1]))
    sex <- vector(length = length(offspring[,1]))
    
    # Analysis 1 : number of heterozygote (he) and homozygote (ho) for each group 
    he5 <- 0 
    ho5 <- 0
    he6 <- 0
    ho6 <- 0
    he7 <- 0
    ho7 <- 0
    
    #Analysis 2 : number of individuals with diplotype identical to father or mother in group 5 and 7 
    m5 <- 0
    f5 <- 0
    m7 <- 0
    fm7 <- 0
    f7 <- 0
    
    # Analysis 3  MHC divergence
    distance <- vector(length = length(offspring[,1]))
    
    # Analysis 4 (haplotype A-H)
    fa <- 0
    fb <- 0
    fc <- 0
    fd <- 0 
    fe <- 0
    ff <- 0
    fg <- 0
    fh <- 0
    
    # Analysis 5 (number of each haplotype inherited from father (mht) or mother (fht))
    fhta <- 0
    fhtb <- 0
    fhtc <- 0
    fhtd <- 0 
    fhte <- 0
    fhtf <- 0
    fhtg <- 0
    fhth <- 0
    mhta <- 0
    mhtb <- 0
    mhtc <- 0
    mhtd <- 0 
    mhte <- 0
    mhtf <- 0
    mhtg <- 0
    mhth <- 0
    
    # Analysis 6 (number of each haplotype received in males (omh) or females (ofh))
    ofha <- 0
    ofhb <- 0
    ofhc <- 0
    ofhd <- 0 
    ofhe <- 0
    ofhf <- 0
    ofhg <- 0
    ofhh <- 0
    omha <- 0
    omhb <- 0
    omhc <- 0
    omhd <- 0 
    omhe <- 0
    omhf <- 0
    omhg <- 0
    omhh <- 0
    
    
    for (i in 1:length(offspring[,1])){
      # haplotype randomly chosen from fathers and mothers
      x <- sample(c(1,2),2,replace=T)
      y <- paste0("f",x[1])
      z <- paste0("m",x[2])
      
      o[i] <- as.character(offspring[i,y])
      p[i] <- as.character(offspring[i,z])
      f1[i] <- as.character(offspring$f1[i])
      f2[i] <- as.character(offspring$f2[i])
      m1[i] <- as.character(offspring$m1[i])
      m2[i] <- as.character(offspring$m2[i])
      sex[i] <- as.character(offspring$sex[i])
      
      # Analysis 3: get MHC divergence of simulated diplotypes
      distance[i] <- div[which(div$h1==o[i]&div$h2==p[i]),3]
      
      # Analysis 1 and 2 
      if (offspring$type[i]==5){
        if (o[i]==p[i]) {
          ho5 <- ho5 + 1    
          if (f1[i]==f2[i]){f5 <- f5 +1}
          else{m5 <- m5 + 1}
          
        } else {
          if (f1[i]!=f2[i]){f5 <-f5 +1}
          else{m5 <- m5 +1 }
          he5 <- he5 + 1 }
      }
      if (offspring$type[i]==6){
        if (o[i]==p[i]) {
          ho6 <- ho6 +1
          
        } else {
          he6 <- he6 + 1 }
      }
      if (offspring$type[i]==7){
        if (o[i]==p[i]) {
          ho7 <- ho7 +1
        } else {
          n <- 0
          m <- 0
          if (f1[i]==o[i]){n <- n +1}
          if (f2[i]==o[i]){n <- n +1}
          if (f1[i]==p[i]){n <- n +1}
          if (f2[i]==p[i]){n <- n +1}
          if (m1[i]==o[i]){m <- m +1}
          if (m2[i]==o[i]){m <- m +1}
          if (m1[i]==p[i]){m <- m +1}
          if (m2[i]==p[i]){m <- m +1}
          fz1[i] <- n
          mz1[i] <- m
          if (n == 1 && m == 1){fm7 <- fm7 + 1}
          if (m == 2){m7 <- m7 + 1}
          if (n == 2){f7 <- f7 + 1}
          he7 <- he7 + 1 }
      }
      
      #Analysis 4 
      if (o[i]=="A") {
        fa <- fa + 1
      } 
      
      if (p[i]=="A") {
        fa <- fa + 1
      }
      if (o[i]=="B") {
        fb <- fb + 1
      } 
      
      if (p[i]=="B") {
        fb <- fb + 1
      }
      if (o[i]=="C") {
        fc <- fc + 1
      } 
      
      if (p[i]=="C") {
        fc <- fc + 1
      }
      if (o[i]=="D") {
        fd <- fd + 1
      } 
      
      if (p[i]=="D") {
        fd <- fd + 1
      }
      if (o[i]=="E") {
        fe <- fe + 1
      } 
      
      if (p[i]=="E") {
        fe <- fe + 1
      }
      if (o[i]=="F") {
        ff <- ff + 1
      } 
      
      if (p[i]=="F") {
        ff <- ff + 1
      }
      if (o[i]=="G") {
        fg <- fg + 1
      } 
      
      if (p[i]=="G") {
        fg <- fg + 1
      }
      if (o[i]=="H") {
        fh <- fh + 1
      } 
      
      if (p[i]=="H") {
        fh <- fh + 1
      }
      
      # Analysis 6 
      if(sex[i] == 1){
        if (o[i]=="A") {
          ofha <- ofha + 1
        } 
        
        if (p[i]=="A") {
          ofha <- ofha + 1
        }
        if (o[i]=="B") {
          ofhb <- ofhb + 1
        } 
        
        if (p[i]=="B") {
          ofhb <- ofhb + 1
        }
        if (o[i]=="C") {
          ofhc <- ofhc + 1
        } 
        
        if (p[i]=="C") {
          ofhc <- ofhc + 1
        }
        if (o[i]=="D") {
          ofhd <- ofhd + 1
        } 
        
        if (p[i]=="D") {
          ofhd <- ofhd + 1
        }
        if (o[i]=="E") {
          ofhe <- ofhe + 1
        } 
        
        if (p[i]=="E") {
          ofhe <- ofhe + 1
        }
        if (o[i]=="F") {
          ofhf <- ofhf + 1
        } 
        
        if (p[i]=="F") {
          ofhf <- ofhf + 1
        }
        if (o[i]=="G") {
          ofhg <- ofhg + 1
        } 
        
        if (p[i]=="G") {
          ofhg <- ofhg + 1
        }
        if (o[i]=="H") {
          ofhh <- ofhh + 1
        } 
        
        if (p[i]=="H") {
          ofhh <- ofhh + 1
        }
        
      }
      if(sex[i] == 2){
        if (o[i]=="A") {
          omha <- omha + 1
        } 
        
        if (p[i]=="A") {
          omha <- omha + 1
        }
        if (o[i]=="B") {
          omhb <- omhb + 1
        } 
        
        if (p[i]=="B") {
          omhb <- omhb + 1
        }
        if (o[i]=="C") {
          omhc <- omhc + 1
        } 
        
        if (p[i]=="C") {
          omhc <- omhc + 1
        }
        if (o[i]=="D") {
          omhd <- omhd + 1
        } 
        
        if (p[i]=="D") {
          omhd <- omhd + 1
        }
        if (o[i]=="E") {
          omhe <- omhe + 1
        } 
        
        if (p[i]=="E") {
          omhe <- omhe + 1
        }
        if (o[i]=="F") {
          omhf <- omhf + 1
        } 
        
        if (p[i]=="F") {
          omhf <- omhf + 1
        }
        if (o[i]=="G") {
          omhg <- omhg + 1
        } 
        
        if (p[i]=="G") {
          omhg <- omhg + 1
        }
        if (o[i]=="H") {
          omhh <- omhh + 1
        } 
        
        if (p[i]=="H") {
          omhh <- omhh + 1
        }
        
      }
      
      # Analysis 5 
      if (offspring$type[i]!= 6){
        if (o[i]=="A") {
          fhta <- fhta + 1
        } 
        
        if (p[i]=="A") {
          mhta <- mhta + 1
        }
        if (o[i]=="B") {
          fhtb <- fhtb + 1
        } 
        
        if (p[i]=="B") {
          mhtb <- mhtb + 1
        }
        if (o[i]=="C") {
          fhtc <- fhtc + 1
        } 
        
        if (p[i]=="C") {
          mhtc <- mhtc + 1
        }
        if (o[i]=="D") {
          fhtd <- fhtd + 1
        } 
        
        if (p[i]=="D") {
          mhtd <- mhtd + 1
        }
        if (o[i]=="E") {
          fhte <- fhte + 1
        } 
        
        if (p[i]=="E") {
          mhte <- mhte + 1
        }
        if (o[i]=="F") {
          fhtf <- fhtf + 1
        } 
        
        if (p[i]=="F") {
          mhtf <- mhtf + 1
        }
        if (o[i]=="G") {
          fhtg <- fhtg + 1
        } 
        
        if (p[i]=="G") {
          mhtg <- mhtg + 1
        }
        if (o[i]=="H") {
          fhth <- fhth + 1
        } 
        
        if (p[i]=="H") {
          mhth <- mhth + 1
        }
      }
    }
    dist[j] <- mean(distance)
    distm[j] <- median(distance)
    he51[j] <- as.numeric(he5)
    ho51[j] <- as.numeric(ho5)
    he61[j] <- as.numeric(he6)
    ho61[j] <- as.numeric(ho6)
    he71[j] <- as.numeric(he7)
    ho71[j] <- as.numeric(ho7)
    m51[j] <- as.numeric(m5)
    f51[j] <- as.numeric(f5)
    m71[j] <- as.numeric(m7)
    f71[j] <- as.numeric(f7)
    fm71[j] <- as.numeric(fm7)
    h5[j] <- as.numeric(he51[j]/ho51[j])
    h6[j] <- as.numeric(he61[j]/ho61[j])
    h7[j] <- as.numeric(he71[j]/ho71[j])
    fga1[j] <- fa
    fgb1[j] <- fb
    fgc1[j] <- fc
    fgd1[j] <- fd
    fge1[j] <- fe
    fgf1[j] <- ff
    fgg1[j] <- fg
    fgh1[j] <- fh
    fhta1[j] <- fhta
    fhtb1[j] <- fhtb
    fhtc1[j] <- fhtc
    fhtd1[j] <- fhtd
    fhte1[j] <- fhte
    fhtf1[j] <- fhtf
    fhtg1[j] <- fhtg
    fhth1[j] <- fhth
    mhta1[j] <- mhta
    mhtb1[j] <- mhtb
    mhtc1[j] <- mhtc
    mhtd1[j] <- mhtd
    mhte1[j] <- mhte
    mhtf1[j] <- mhtf
    mhtg1[j] <- mhtg
    mhth1[j] <- mhth
    ofha1[j] <- ofha
    ofhb1[j] <- ofhb
    ofhc1[j] <- ofhc
    ofhd1[j] <- ofhd
    ofhe1[j] <- ofhe
    ofhf1[j] <- ofhf
    ofhg1[j] <- ofhg
    ofhh1[j] <- ofhh
    omha1[j] <- omha
    omhb1[j] <- omhb
    omhc1[j] <- omhc
    omhd1[j] <- omhd
    omhe1[j] <- omhe
    omhf1[j] <- omhf
    omhg1[j] <- omhg
    omhh1[j] <- omhh
  }
  mylist <- list("dist"=dist,"m51"=m51,"f51"=f51,"m71"=m71,"f71"=f71,"fm71"=fm71,"h5"=h5,"h6"=h6,"h7"=h7, "fga1"=fga1,"fgb1"=fgb1,"fgc1"=fgc1,"fgd1"=fgd1,"fge1"=fge1,"fgf1"=fgf1,"fgg1"=fgg1,"fgh1"=fgh1,"fhta1"=fhta1,"fhtb1"=fhtb1,"fhtc1"=fhtc1,"fhtd1"=fhtd1,"fhte1"=fhte1,"fhtf1"=fhtf1,"fhtg1"=fhtg1,"fhth1"=fhth1,"mhta1"=mhta1,"mhtb1"=mhtb1,"mhtc1"=mhtc1,"mhtd1"=mhtd1,"mhte1"=mhte1,"mhtf1"=mhtf1,"mhtg1"=mhtg1,"mhth1"=mhth1,"ofha1"=ofha1,"ofhb1"=ofhb1,"ofhc1"=ofhc1,"ofhd1"=ofhd1,"ofhe1"=ofhe1,"ofhf1"=ofhf1,"ofhg1"=ofhg1,"ofhh1"=ofhh1,"omha1"=omha1,"omhb1"=omhb1,"omhc1"=omhc1,"omhd1"=omhd1,"omhe1"=omhe1,"omhf1"=omhf1,"omhg1"=omhg1,"omhh1"=omhh1)
  return(mylist)
}

#choose the number of simualation iteration 
result <- simulation_post(10000)

#Real data
#haplotype of offspring(o/p), mother(f1/f2) and father(m1/m2)
o <- vector(length = length(offspring[,1]))
p <- vector(length = length(offspring[,1]))
f1 <- vector(length = length(offspring[,1]))
m1 <- vector(length = length(offspring[,1]))
f2 <- vector(length = length(offspring[,1]))
m2 <- vector(length = length(offspring[,1]))

#Analysis 5 haplotype inherited from mother (rfht) and father(rmht)
rfht <- vector(length = length(offspring[,1]))
rmht <- vector(length = length(offspring[,1]))

#Analysis 3 :MHC divergence
rdistance <- vector(length = length(offspring[,1]))
rdist <- 0

#Analysis 1
rhe1 <- 0 
rho1 <- 0
rhe2 <- 0
rho2 <- 0
rhe3 <- 0 
rho3 <- 0
rhe4 <- 0
rho4 <- 0
rhe5 <- 0 
rho5 <- 0
rhe6 <- 0
rho6 <- 0
rhe7 <- 0
rho7 <- 0

#Analysis 2
rm5 <- 0
rf5 <- 0
rm7 <- 0
rfm7 <- 0
rf7 <- 0

#Analysis 4
rfa <- 0
rfb <- 0
rfc <- 0
rfd <- 0 
rfe <- 0
rff <- 0
rfg <- 0
rfh <- 0

#Analysis 6 : MHC haplotype in females (rof) and males (rom)
rofa <- 0
rofb <- 0
rofc <- 0
rofd <- 0 
rofe <- 0
roff <- 0
rofg <- 0
rofh <- 0
roma <- 0
romb <- 0
romc <- 0
romd <- 0 
rome <- 0
romf <- 0
romg <- 0
romh <- 0

for (i in 1:length(offspring[,1])){
  
  o[i] <- as.character(offspring[i,"h1"])
  p[i] <- as.character(offspring[i,"h2"])
  f1[i] <- as.character(offspring[i,"f1"])
  f2[i] <- as.character(offspring[i,"f2"])
  m1[i] <- as.character(offspring[i,"m1"])
  m2[i] <- as.character(offspring[i,"m2"])
  rdistance[i] <- div[which(div$h1==o[i]&div$h2==p[i]),3]
  if (offspring$type[i]==3){
    if (f1[i]==f2[i]&f1[i]==o[i]) {
      rfht[i] <- o[i]
      rmht[i] <- p[i]
    } 
    if (f1[i]==f2[i]&f1[i]==p[i]) {
      rmht[i] <- o[i]
      rfht[i] <- p[i]
    }
    if (m1[i]==m2[i]&m1[i]==p[i]) {
      rfht[i] <- o[i]
      rmht[i] <- p[i]
    } 
    if (m1[i]==m2[i]&m1[i]==o[i]) {
      rfht[i] <- p[i]
      rmht[i] <- o[i]
    } 
  }
  if (offspring$type[i]==4){
    if (o[i]==f1[i]|o[i]==f2[i]) {
      rfht[i] <- o[i]
      rmht[i] <- p[i]
    } else {
      rmht[i] <- o[i]
      rfht[i] <- p[i]
    }
  }
  if (offspring$type[i]==5){
    if (o[i]==p[i]) {
      rfht[i] <- o[i]
      rmht[i] <- p[i]
      rho5 <- rho5 + 1    
      if (f1[i]==f2[i]){rf5 <- rf5 +1}
      else{rm5 <- rm5 + 1}
      
    } else {
      if (f1[i]==f2[i]){
        rm5 <-rm5 +1
        if(o[i]==f1[i]){
          rfht[i] <- o[i]
          rmht[i] <- p[i]
        }else{
          rmht[i] <- o[i]
          rfht[i] <- p[i]}
      }else{
        rf5 <- rf5 +1
        if(p[i]==m1[i]){
          rfht[i] <- o[i]
          rmht[i] <- p[i]
        }else{
          rmht[i] <- o[i]
          rfht[i] <- p[i]}
      }
      rhe5 <- rhe5 + 1 }
  }
  if (offspring$type[i]==6){
    if (o[i]==p[i]) {
      rho6 <- rho6 +1
      
    } else {
      rhe6 <- rhe6 + 1 }
  }
  if (offspring$type[i]==7){
    if (o[i]==p[i]) {
      rho7 <- rho7 +1
      rfht[i] <- o[i]
      rmht[i] <- p[i]
    } else {
      n <- 0
      m <- 0
      if (f1[i]==o[i]){n <- n +1}
      if (f2[i]==o[i]){n <- n +1}
      if (f1[i]==p[i]){n <- n +1}
      if (f2[i]==p[i]){n <- n +1}
      if (m1[i]==o[i]){m <- m +1}
      if (m2[i]==o[i]){m <- m +1}
      if (m1[i]==p[i]){m <- m +1}
      if (m2[i]==p[i]){m <- m +1}
      if (n ==1 && m == 1){rfm7 <- rfm7 + 1
      if (o[i]==f1[i]|o[i]==f2[i]) {
        rfht[i] <- o[i]
        rmht[i] <- p[i]
      } else {
        rmht[i] <- o[i]
        rfht[i] <- p[i]}
      }
      
      if (m ==2){rm7 <- rm7 + 1
      if (o[i]!=m1[i]|o[i]!=m2[i]) {
        rfht[i] <- o[i]
        rmht[i] <- p[i]
      } else {
        rmht[i] <- o[i]
        rfht[i] <- p[i]}
      }
      
      if (n == 2){rf7 <- rf7 + 1
      if (o[i]!=f1[i]|o[i]!=f2[i]) {
        rfht[i] <- o[i]
        rmht[i] <- p[i]
      } else {
        rmht[i] <- o[i]
        rfht[i] <- p[i]}
      }
      rhe7 <- rhe7 + 1 }
  }
  
  if (o[i]=="A") {
    rfa <- rfa + 1
  } 
  
  if (p[i]=="A") {
    rfa <- rfa + 1
  }
  if (o[i]=="B") {
    rfb <- rfb + 1
  } 
  
  if (p[i]=="B") {
    rfb <- rfb + 1
  }
  if (o[i]=="C") {
    rfc <- rfc + 1
  } 
  
  if (p[i]=="C") {
    rfc <- rfc + 1
  }
  if (o[i]=="D") {
    rfd <- rfd + 1
  } 
  
  if (p[i]=="D") {
    rfd <- rfd + 1
  }
  if (o[i]=="E") {
    rfe <- rfe + 1
  } 
  
  if (p[i]=="E") {
    rfe <- rfe + 1
  }
  if (o[i]=="F") {
    rff <- rff + 1
  } 
  
  if (p[i]=="F") {
    rff <- rff + 1
  }
  if (o[i]=="G") {
    rfg <- rfg + 1
  } 
  
  if (p[i]=="G") {
    rfg <- rfg + 1
  }
  if (o[i]=="H") {
    rfh <- rfh + 1
  } 
  
  if (p[i]=="H") {
    rfh <- rfh + 1
  }
  
  if (offspring$sex[i]==1){
    if (o[i]=="A") {
      rofa <- rofa + 1
    } 
    
    if (p[i]=="A") {
      rofa <- rofa + 1
    }
    if (o[i]=="B") {
      rofb <- rofb + 1
    } 
    
    if (p[i]=="B") {
      rofb <- rofb + 1
    }
    if (o[i]=="C") {
      rofc <- rofc + 1
    } 
    
    if (p[i]=="C") {
      rofc <- rofc + 1
    }
    if (o[i]=="D") {
      rofd <- rofd + 1
    } 
    
    if (p[i]=="D") {
      rofd <- rofd + 1
    }
    if (o[i]=="E") {
      rofe <- rofe + 1
    } 
    
    if (p[i]=="E") {
      rofe <- rofe + 1
    }
    if (o[i]=="F") {
      roff <- roff + 1
    } 
    
    if (p[i]=="F") {
      roff <- roff + 1
    }
    if (o[i]=="G") {
      rofg <- rofg + 1
    } 
    
    if (p[i]=="G") {
      rofg <- rofg + 1
    }
    if (o[i]=="H") {
      rofh <- rofh + 1
    } 
    
    if (p[i]=="H") {
      rofh <- rofh + 1
    }
  }
  
  if (offspring$sex[i]==2){
    if (o[i]=="A") {
      roma <- roma + 1
    } 
    
    if (p[i]=="A") {
      roma <- roma + 1
    }
    if (o[i]=="B") {
      romb <- romb + 1
    } 
    
    if (p[i]=="B") {
      romb <- romb + 1
    }
    if (o[i]=="C") {
      romc <- romc + 1
    } 
    
    if (p[i]=="C") {
      romc <- romc + 1
    }
    if (o[i]=="D") {
      romd <- romd + 1
    } 
    
    if (p[i]=="D") {
      romd <- romd + 1
    }
    if (o[i]=="E") {
      rome <- rome + 1
    } 
    
    if (p[i]=="E") {
      rome <- rome + 1
    }
    if (o[i]=="F") {
      romf <- romf + 1
    } 
    
    if (p[i]=="F") {
      romf <- romf + 1
    }
    if (o[i]=="G") {
      romg <- romg + 1
    } 
    
    if (p[i]=="G") {
      romg <- romg + 1
    }
    if (o[i]=="H") {
      romh <- romh + 1
    } 
    
    if (p[i]=="H") {
      romh <- romh + 1
    }
  }
}
rh5 <- as.numeric(rhe5/rho5)
rh6 <- as.numeric(rhe6/rho6)
rh7 <- as.numeric(rhe7/rho7)
rdist <- mean(rdistance)

offspring$rfht <- rfht
offspring$rmht <- rmht
offspring1 <- offspring[offspring$type!=6,]
o <- vector(length = length(offspring1[,1]))
p <- vector(length = length(offspring1[,1]))
rfhta <- 0
rfhtb <- 0
rfhtc <- 0
rfhtd <- 0
rfhte <- 0
rfhtf <- 0
rfhtg <- 0
rfhth <- 0
rmhta <- 0
rmhtb <- 0
rmhtc <- 0
rmhtd <- 0
rmhte <- 0
rmhtf <- 0
rmhtg <- 0
rmhth <- 0


for (i in 1:length(offspring1[,1])){
  o[i] <- offspring1$rfht[i]
  p[i] <- offspring1$rmht[i]
  if (o[i]=="A") {
    rfhta <- rfhta + 1
  } 
  
  if (p[i]=="A") {
    rmhta <- rmhta + 1
  }
  if (o[i]=="B") {
    rfhtb <- rfhtb + 1
  } 
  
  if (p[i]=="B") {
    rmhtb <- rmhtb + 1
  }
  if (o[i]=="C") {
    rfhtc <- rfhtc + 1
  } 
  
  if (p[i]=="C") {
    rmhtc <- rmhtc + 1
  }
  if (o[i]=="D") {
    rfhtd <- rfhtd + 1
  } 
  
  if (p[i]=="D") {
    rmhtd <- rmhtd + 1
  }
  if (o[i]=="E") {
    rfhte <- rfhte + 1
  } 
  
  if (p[i]=="E") {
    rmhte <- rmhte + 1
  }
  if (o[i]=="F") {
    rfhtf <- rfhtf + 1
  } 
  
  if (p[i]=="F") {
    rmhtf <- rmhtf + 1
  }
  if (o[i]=="G") {
    rfhtg <- rfhtg + 1
  } 
  
  if (p[i]=="G") {
    rmhtg <- rmhtg + 1
  }
  if (o[i]=="H") {
    rfhth <- rfhth + 1
  } 
  
  if (p[i]=="H") {
    rmhth <- rmhth + 1
  }
}




# total number of MHC haplotypes 
fg1 <- sum(rfa,rfb,rfc,rfd,rfe,rff,rfg,rfh)
fht1 <- sum(rfhta,rfhtb,rfhtc,rfhtd,rfhte,rfhtf,rfhtg,rfhth)
mht1 <- sum(rmhta,rmhtb,rmhtc,rmhtd,rmhte,rmhtf,rmhtg,rmhth)
ofg1 <- sum(rofa,rofb,rofc,rofd,rofe,roff,rofg,rofh)
omg1 <- sum(roma,romb,romc,romd,rome,romf,romg,romh)


h1 <- result$h1
h2 <- result$h2
h3 <- result$h3
h4 <- result$h4
h5 <- result$h5
h6 <- result$h6
h7 <- result$h7

he11 <- result$he11
ho11 <- result$ho11
he21 <- result$he21
ho21 <- result$ho21
he31 <- result$he31
ho31 <- result$ho31
he41 <- result$he41
ho41 <- result$ho41
he51 <- result$he51
ho51 <- result$ho51
he61 <- result$he51
ho61 <- result$ho61
he71 <- result$he71
ho71 <- result$ho71
f51 <- result$f51
m51 <- result$m51
f71 <- result$f71
m71 <- result$m71
fm71 <- result$fm71
fga1 <- result$fga1
fgb1 <- result$fgb1
fgc1 <- result$fgc1
fgd1 <- result$fgd1
fge1 <- result$fge1
fgf1 <- result$fgf1
fgg1 <- result$fgg1
fgh1 <- result$fgh1
fhta1 <- result$fhta1
fhtb1 <- result$fhtb1
fhtc1 <- result$fhtc1
fhtd1 <- result$fhtd1
fhte1 <- result$fhte1
fhtf1 <- result$fhtf1
fhtg1 <- result$fhtg1
fhth1 <- result$fhth1
mhta1 <- result$mhta1
mhtb1 <- result$mhtb1
mhtc1 <- result$mhtc1
mhtd1 <- result$mhtd1
mhte1 <- result$mhte1
mhtf1 <- result$mhtf1
mhtg1 <- result$mhtg1
mhth1 <- result$mhth1
ofha1 <- result$ofha1
ofhb1 <- result$oghb1
ofhc1 <- result$ofhc1
ofhd1 <- result$ofhd1
ofhe1 <- result$ofhe1
ofhf1 <- result$ofhf1
ofhg1 <- result$ofhg1
ofhh1 <- result$ofhh1
omha1 <- result$omha1
omhb1 <- result$omhb1
omhc1 <- result$omhc1
omhd1 <- result$omhd1
omhe1 <- result$omhe1
omhf1 <- result$omhf1
omhg1 <- result$omhg1
omhh1 <- result$omhh1

he11 <- result$he11
ho11 <- result$ho11
he21 <- result$he21
ho21 <- result$ho21
he31 <- result$he31
ho31 <- result$ho31
he41 <- result$he41
ho41 <- result$ho41
he51 <- result$he51
ho51 <- result$ho51
he61 <- result$he51
ho61 <- result$ho61
he71 <- result$he71
ho71 <- result$ho71
f51 <- result$f51
m51 <- result$m51
f71 <- result$f71
m71 <- result$m71
fm71 <- result$fm71
fga1 <- result$fga1
fgb1 <- result$fgb1
fgc1 <- result$fgc1
fgd1 <- result$fgd1
fge1 <- result$fge1
fgf1 <- result$fgf1
fgg1 <- result$fgg1
fgh1 <- result$fgh1
fhta1 <- result$fhta1
fhtb1 <- result$fhtb1
fhtc1 <- result$fhtc1
fhtd1 <- result$fhtd1
fhte1 <- result$fhte1
fhtf1 <- result$fhtf1
fhtg1 <- result$fhtg1
fhth1 <- result$fhth1
mhta1 <- result$mhta1
mhtb1 <- result$mhtb1
mhtc1 <- result$mhtc1
mhtd1 <- result$mhtd1
mhte1 <- result$mhte1
mhtf1 <- result$mhtf1
mhtg1 <- result$mhtg1
mhth1 <- result$mhth1
ofha1 <- result$ofha1
ofhb1 <- result$ofhb1
ofhc1 <- result$ofhc1
ofhd1 <- result$ofhd1
ofhe1 <- result$ofhe1
ofhf1 <- result$ofhf1
ofhg1 <- result$ofhg1
ofhh1 <- result$ofhh1
omha1 <- result$omha1
omhb1 <- result$omhb1
omhc1 <- result$omhc1
omhd1 <- result$omhd1
omhe1 <- result$omhe1
omhf1 <- result$omhf1
omhg1 <- result$omhg1
omhh1 <- result$omhh1
dist <-result$dist

# make figures
data1 <- data.frame(h5,h6,h7)
data2 <- data.frame(f51,m51,f71,m71,fm71)
data3 <- data.frame(fga1,fgb1,fgc1,fgd1,fge1,fgf1,fgg1,fgh1)
data4 <- data.frame(fhta1,fhtb1,fhtc1,fhtd1,fhte1,fhtf1,fhtg1,fhth1)
data5 <- data.frame(mhta1,mhtb1,mhtc1,mhtd1,mhte1,mhtf1,mhtg1,mhth1)
data6 <- data.frame(ofha1,ofhb1,ofhc1,ofhd1,ofhe1,ofhf1,ofhg1,ofhh1)
data7 <- data.frame(omha1,omhb1,omhc1,omhd1,omhe1,omhf1,omhg1,omhh1)
data8 <- data.frame(dist)

#analysis 1
p1 <- ggplot(aes(x=h5), data=data1) + geom_vline(xintercept=rh5,lty=2,color="blue")+geom_histogram(binwidth = 0.1,alpha=0.6) + scale_x_continuous(breaks = seq(0, 4, by = 0.5))+xlab("Ratio of heterozygote/homozygote")+ylab("Number of simulations")+ggtitle("Group 5")+geom_vline(xintercept=quantile(h5,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(h5,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())
p2 <- ggplot(aes(x=h6), data=data1) + geom_vline(xintercept=rh6,lty=2,color="blue")+geom_histogram(binwidth = 0.1,alpha=0.6) + scale_x_continuous(breaks = seq(0, 4, by = 0.5))+xlab("Ratio of heterozygote/homozygote")+ylab("Number of simulations")+ggtitle("Group 6")+geom_vline(xintercept=quantile(h6,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(h6,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())
p3 <- ggplot(aes(x=h7), data=data1) + geom_vline(xintercept=rh7,lty=2,color="blue")+geom_histogram(binwidth = 0.1,alpha=0.6) + scale_x_continuous(breaks = seq(0, 4, by = 0.5))+xlab("Ratio of heterozygote/homozygote")+ylab("Number of simulations")+ggtitle("Group 7")+geom_vline(xintercept=quantile(h7,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(h7,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())

#analysis 2
p4 <- ggplot(aes(x=f51), data=data2) + geom_vline(xintercept=rf5,lty=2,color="blue")+geom_histogram(binwidth = 1,alpha=0.6) + scale_x_continuous(breaks = seq(0, 2000, by =10))+xlab("Number of individuals with maternal identical diplotype")+ylab("Number of simulations")+ggtitle("Group 5")+geom_vline(xintercept=quantile(f51,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(f51,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())
p5 <- ggplot(aes(x=m51), data=data2) + geom_vline(xintercept=rm5,lty=2,color="blue")+geom_histogram(binwidth = 1,alpha=0.6) + scale_x_continuous(breaks = seq(0, 2000, by = 10))+xlab("Number of individuals with paternal identical diplotype")+ylab("Number of simulations")+ggtitle("Group 5")+geom_vline(xintercept=quantile(m51,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(m51,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())
p6 <- ggplot(aes(x=f71), data=data2) + geom_vline(xintercept=rf7,lty=2,color="blue")+geom_histogram(binwidth = 1,alpha=0.6) + scale_x_continuous(breaks = seq(0, 2000, by = 20))+xlab("Number of individuals with maternal identical diplotype")+ylab("Number of simulations")+ggtitle("Group 7")+geom_vline(xintercept=quantile(f71,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(f71,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())
ggplot(aes(x=fm71), data=data2) + geom_vline(xintercept=rfm7,lty=2,color="blue")+geom_histogram(binwidth = 1,alpha=0.6) + scale_x_continuous(breaks = seq(0, 2000, by = 20))+geom_vline(xintercept=quantile(fm71,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(fm71,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())
p7 <- ggplot(aes(x=m71), data=data2) + geom_vline(xintercept=rm7,lty=2,color="blue")+geom_histogram(binwidth = 1,alpha=0.6) + scale_x_continuous(breaks = seq(0, 2000, by = 20))+xlab("Number of individuals with paternal identical diplotype")+ylab("Number of simulations")+ggtitle("Group 7")+geom_vline(xintercept=quantile(m71,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(m71,0.025),lty=3,size=1)+theme(axis.title.y = element_blank())

#analysis 3
p8 <- ggplot(aes(x=dist),data=data8) + geom_vline(xintercept=rdist,lty=2,color="blue")+geom_histogram(binwidth = 0.0001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 0.2, by = 0.001))+xlab("Mean p-distance of MHC diplotype across all individuals"  )+ylab("Number of simulations")+geom_vline(xintercept=quantile(dist,0.975),lty=3,size=1)+geom_vline(xintercept=quantile(dist,0.025),lty=3,size=1)

#analysis 4
fa <- ggplot(aes(x=fga1/fg1), data=data3) + geom_vline(xintercept=rfa/fg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype A")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fga1/fg1,0.003125),lty=3,size=1)+geom_vline(xintercept=quantile(fga1/fg1,0.996875),lty=3,size=1)
fb <- ggplot(aes(x=fgb1/fg1), data=data3) + geom_vline(xintercept=rfb/fg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype B")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fgb1/fg1,0.003125),lty=3,size=1)+geom_vline(xintercept=quantile(fgb1/fg1,0.996875),lty=3,size=1)
fc <- ggplot(aes(x=fgc1/fg1), data=data3) + geom_vline(xintercept=rfc/fg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype C")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fgc1/fg1,0.003125),lty=3,size=1)+geom_vline(xintercept=quantile(fgc1/fg1,0.996875),lty=3,size=1)
fd <- ggplot(aes(x=fgd1/fg1), data=data3) + geom_vline(xintercept=rfd/fg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype D")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fgd1/fg1,0.003125),lty=3,size=1)+geom_vline(xintercept=quantile(fgd1/fg1,0.996875),lty=3,size=1)                                                                                       
fe <- ggplot(aes(x=fge1/fg1), data=data3) + geom_vline(xintercept=rfe/fg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype E")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fge1/fg1,0.003125),lty=3,size=1)+geom_vline(xintercept=quantile(fge1/fg1,0.996875),lty=3,size=1)
ff <- ggplot(aes(x=fgf1/fg1), data=data3) + geom_vline(xintercept=rff/fg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype F")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fgf1/fg1,0.003125),lty=3,size=1)+geom_vline(xintercept=quantile(fgf1/fg1,0.996875),lty=3,size=1)
fg <- ggplot(aes(x=fgg1/fg1), data=data3) + geom_vline(xintercept=rfg/fg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype G")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fgg1/fg1,0.003125),lty=3,size=1)+geom_vline(xintercept=quantile(fgg1/fg1,0.996875),lty=3,size=1)
fh <- ggplot(aes(x=fgh1/fg1), data=data3) + geom_vline(xintercept=rfh/fg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype H")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fgh1/fg1,0.003125),lty=3,size=1)+geom_vline(xintercept=quantile(fgh1/fg1,0.996875),lty=3,size=1)


#analysis 5 (Frequency)
fhta <- ggplot(aes(x=fhta1/fht1), data=data4) + geom_vline(xintercept=rfhta/fht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of maternal haplotype A")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fhta1/fht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fhta1/fht1,0.9984375),lty=3,size=1)
fhtb <- ggplot(aes(x=fhtb1/fht1), data=data4) + geom_vline(xintercept=rfhtb/fht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of maternal haplotype B")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fhtb1/fht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fhtb1/fht1,0.9984375),lty=3,size=1)
fhtc <- ggplot(aes(x=fhtc1/fht1), data=data4) + geom_vline(xintercept=rfhtc/fht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of maternal haplotype C")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fhtc1/fht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fhtc1/fht1,0.9984375),lty=3,size=1)
fhtd <- ggplot(aes(x=fhtd1/fht1), data=data4) + geom_vline(xintercept=rfhtd/fht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of maternal haplotype D")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fhtd1/fht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fhtd1/fht1,0.9984375),lty=3,size=1)
fhte <- ggplot(aes(x=fhte1/fht1), data=data4) + geom_vline(xintercept=rfhte/fht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of maternal haplotype E")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fhte1/fht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fhte1/fht1,0.9984375),lty=3,size=1)
fhtf <- ggplot(aes(x=fhtf1/fht1), data=data4) + geom_vline(xintercept=rfhtf/fht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of maternal haplotype F")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fhtf1/fht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fhtf1/fht1,0.9984375),lty=3,size=1)
fhtg <- ggplot(aes(x=fhtg1/fht1), data=data4) + geom_vline(xintercept=rfhtg/fht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of maternal haplotype G")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fhtg1/fht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fhtg1/fht1,0.9984375),lty=3,size=1)
fhth <- ggplot(aes(x=fhth1/fht1), data=data4) + geom_vline(xintercept=rfhth/fht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of maternal haplotype H")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(fhth1/fht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(fhth1/fht1,0.9984375),lty=3,size=1)


mhta <- ggplot(aes(x=mhta1/mht1), data=data5) + geom_vline(xintercept=rmhta/mht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of paternal haplotype A")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(mhta1/mht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mhta1/mht1,0.9984375),lty=3,size=1)
mhtb <- ggplot(aes(x=mhtb1/mht1), data=data5) + geom_vline(xintercept=rmhtb/mht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of paternal haplotype B")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(mhtb1/mht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mhtb1/mht1,0.9984375),lty=3,size=1)
mhtc <- ggplot(aes(x=mhtc1/mht1), data=data5) + geom_vline(xintercept=rmhtc/mht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of paternal haplotype C")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(mhtc1/mht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mhtc1/mht1,0.9984375),lty=3,size=1)
mhtd <- ggplot(aes(x=mhtd1/mht1), data=data5) + geom_vline(xintercept=rmhtd/mht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of paternal haplotype D")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(mhtd1/mht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mhtd1/mht1,0.9984375),lty=3,size=1)
mhte <- ggplot(aes(x=mhte1/mht1), data=data5) + geom_vline(xintercept=rmhte/mht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of paternal haplotype E")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(mhte1/mht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mhte1/mht1,0.9984375),lty=3,size=1)
mhtf <- ggplot(aes(x=mhtf1/mht1), data=data5) + geom_vline(xintercept=rmhtf/mht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of paternal haplotype F")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(mhtf1/mht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mhtf1/mht1,0.9984375),lty=3,size=1)
mhtg <- ggplot(aes(x=mhtg1/mht1), data=data5) + geom_vline(xintercept=rmhtg/mht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of paternal haplotype G")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(mhtg1/mht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mhtg1/mht1,0.9984375),lty=3,size=1)
mhth <- ggplot(aes(x=mhth1/mht1), data=data5) + geom_vline(xintercept=rmhth/mht1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of paternal haplotype H")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(mhth1/mht1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(mhth1/mht1,0.9984375),lty=3,size=1)

#analysis 6 (Frequency)
ofa <- ggplot(aes(x=ofha1/ofg1), data=data6) + geom_vline(xintercept=rofa/ofg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype A in females")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(ofha1/ofg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(ofha1/ofg1,0.9984375),lty=3,size=1)
ofb <- ggplot(aes(x=ofhb1/ofg1), data=data6) + geom_vline(xintercept=rofb/ofg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype B in females")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(ofhb1/ofg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(ofhb1/ofg1,0.9984375),lty=3,size=1)
ofc <- ggplot(aes(x=ofhc1/ofg1), data=data6) + geom_vline(xintercept=rofc/ofg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype C in females")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(ofhc1/ofg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(ofhc1/ofg1,0.9984375),lty=3,size=1)
ofd <- ggplot(aes(x=ofhd1/ofg1), data=data6) + geom_vline(xintercept=rofd/ofg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype D in females")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(ofhd1/ofg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(ofhd1/ofg1,0.9984375),lty=3,size=1)
ofe <- ggplot(aes(x=ofhe1/ofg1), data=data6) + geom_vline(xintercept=rofe/ofg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype E in females")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(ofhe1/ofg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(ofhe1/ofg1,0.9984375),lty=3,size=1)
off <- ggplot(aes(x=ofhf1/ofg1), data=data6) + geom_vline(xintercept=roff/ofg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype F in females")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(ofhf1/ofg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(ofhf1/ofg1,0.9984375),lty=3,size=1)
ofg <- ggplot(aes(x=ofhg1/ofg1), data=data6) + geom_vline(xintercept=rofg/ofg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype G in females")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(ofhg1/ofg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(ofhg1/ofg1,0.9984375),lty=3,size=1)
ofh <- ggplot(aes(x=ofhh1/ofg1), data=data6) + geom_vline(xintercept=rofh/ofg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype H in females")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(ofhh1/ofg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(ofhh1/ofg1,0.9984375),lty=3,size=1)

oma <- ggplot(aes(x=omha1/omg1), data=data7) + geom_vline(xintercept=roma/omg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype A in males")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(omha1/omg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(omha1/omg1,0.9984375),lty=3,size=1)
omb <- ggplot(aes(x=omhb1/omg1), data=data7) + geom_vline(xintercept=romb/omg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype B in males")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(omhb1/omg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(omhb1/omg1,0.9984375),lty=3,size=1)
omc <- ggplot(aes(x=omhc1/omg1), data=data7) + geom_vline(xintercept=romc/omg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype C in males")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(omhc1/omg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(omhc1/omg1,0.9984375),lty=3,size=1)
omd <- ggplot(aes(x=omhd1/omg1), data=data7) + geom_vline(xintercept=romd/omg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype D in males")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(omhd1/omg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(omhd1/omg1,0.9984375),lty=3,size=1)
ome <- ggplot(aes(x=omhe1/omg1), data=data7) + geom_vline(xintercept=rome/omg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype E in males")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(omhe1/omg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(omhe1/omg1,0.9984375),lty=3,size=1)
omf <- ggplot(aes(x=omhf1/omg1), data=data7) + geom_vline(xintercept=romf/omg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype F in males")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(omhf1/omg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(omhf1/omg1,0.9984375),lty=3,size=1)
omg <- ggplot(aes(x=omhg1/omg1), data=data7) + geom_vline(xintercept=romg/omg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype G in males")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(omhg1/omg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(omhg1/omg1,0.9984375),lty=3,size=1)
omh <- ggplot(aes(x=omhh1/omg1), data=data7) + geom_vline(xintercept=romh/omg1,lty=2,color="blue")+geom_histogram(binwidth = 0.001,alpha=0.6) + scale_x_continuous(breaks = seq(0, 1, by = 0.01))+xlab("Frequency of haplotype H in males")+theme(axis.title.y = element_blank())+geom_vline(xintercept=quantile(omhh1/omg1,0.0015625),lty=3,size=1)+geom_vline(xintercept=quantile(omhh1/omg1,0.9984375),lty=3,size=1)


fg31 <- ggarrange(p1,p2,p3,nrow=1,ncol=3)
fg41 <- ggarrange(p4,p5,p6,p7,nrow=2,ncol=2)

fgeno <- ggarrange(fa,fb,fc,fd,fe,ff,fg,fh,nrow=4,ncol=2)
fht <- ggarrange(fhta,fhtb,fhtc,fhtd,fhte,fhtf,fhtg,fhth,nrow=4,ncol=2)
mht <- ggarrange(mhta,mhtb,mhtc,mhtd,mhte,mhtf,mhtg,mhth,nrow=4,ncol=2)
of <- ggarrange(ofa,ofb,ofc,ofd,ofe,off,ofg,ofh,nrow=4,ncol=2)
om <- ggarrange(oma,omb,omc,omd,ome,omf,omg,omh,nrow=4,ncol=2)

annotate_figure(fgeno,left = text_grob("Number of simulations", rot = 90))
annotate_figure(fg31,left = text_grob("Number of simulations", rot = 90))
annotate_figure(fg41,left = text_grob("Number of simulations", rot = 90))


fht <- annotate_figure(fht,left = text_grob("Number of simulations", rot = 90))
mht <- annotate_figure(mht,left = text_grob("Number of simulations", rot = 90))
of <- annotate_figure(of,left = text_grob("Number of simulations", rot = 90))
om <- annotate_figure(om,left = text_grob("Number of simulations", rot = 90))
ggarrange(fht,of,labels=c("A","B"))
ggarrange(mht,om,labels=c("A","B"))

