Define parameters, functions:

numinds=1000
numloci=100
kneecap<-function(x){ if(x>1){ return(1) } else(return(x))} # turn loci from additive > dominant (two 1's == one 1)
library(ggplot2)
hetplot=function(a,b){  
  bob=data.frame(a,b)
  namex=deparse(substitute(a));
  namey=deparse(substitute(b));
  the_plot=ggplot(bob)+geom_point(aes(y=b,x=a)) + 
    geom_smooth(aes(y=b,x=a),method="lm") +   
    geom_text(aes(-Inf,Inf,label=paste("r2=",round(summary(lm(b~a))$r.squared,2))),vjust=2,hjust=-0.5) +
    xlab(namex)+ylab(namey);
 print(the_plot)
}

Make parents, F1

Inbreds and phenotypes

There’s no dominance in inbreds, so we just sum to get phenotypes.

p1 <- lapply(1:numinds, function(X) sample(c(0,1),numloci,replace=T))
p1_phenotype=unlist(lapply(p1,sum))
p2 <- lapply(1:numinds, function(X) sample(c(0,1),numloci,replace=T))
p2_phenotype=unlist(lapply(p2,sum))

Identify best parent

best_p<-sapply(1:1000, function(i) max(p1_phenotype[i],p2_phenotype[i]))

Make F1

f1 <- lapply(1:numinds, function(X) p1[[X]]+p2[[X]])

Make F1 phenotypes assuming additiviy, dominance, or both

# purely additive phenotype for the F1
f1_a_phenotype<- unlist(lapply(f1,sum))

#purely dominant phenotype for the F1
d1=f1
d1<-lapply(1:numinds, function(X) sapply(1:numloci, function(Y) kneecap(d1[[X]][Y]) ))
f1_d_phenotype<-unlist(lapply(d1,sum))

#both, with some loci dominant
domloci<-sample(1:numloci,round(runif(1,1,numloci))) # random sample of loci that we will make dominant; rest are additive
print(length(domloci))
## [1] 27
ad1=f1
ad1<-lapply(1:numinds, function(X) sapply(1:numloci, function(Y) if(Y %in% domloci){ kneecap(ad1[[X]][Y]) } else{ return(ad1[[X]][Y])} ))
f1_ad_phenotype<-unlist(lapply(ad1,sum))

Calculate heterosis for all three F1 types

# heterosis additive
HPH_a<-f1_a_phenotype-best_p
pHPH_a=HPH_a/best_p

#heterosis dominant
HPH_d<-f1_d_phenotype-best_p
pHPH_d<-HPH_d/best_p

#heterosis both
HPH_ad<-f1_ad_phenotype-best_p
pHPH_ad<-HPH_ad/best_p

Plot correlations

Heterosis is caused by additive and dominance.

We map as additive:

hetplot(f1_a_phenotype,HPH_a)

hetplot(f1_a_phenotype,pHPH_a)

We map as both:

hetplot(f1_ad_phenotype,HPH_a)

hetplot(f1_ad_phenotype,pHPH_a)

We map as dominance:

hetplot(f1_d_phenotype,HPH_ad)

hetplot(f1_d_phenotype,pHPH_ad)

Heterosis is caused by dominance.

We map as additive:

hetplot(f1_a_phenotype,HPH_d)

hetplot(f1_a_phenotype,pHPH_d)

We map as both:

hetplot(f1_ad_phenotype,HPH_d)

hetplot(f1_ad_phenotype,pHPH_d)

We map as dominance:

hetplot(f1_d_phenotype,HPH_d)

hetplot(f1_d_phenotype,pHPH_d)

Heterosis is caused by additive and dominance.

We map as additive:

hetplot(f1_a_phenotype,HPH_ad)

hetplot(f1_a_phenotype,pHPH_ad)

We map as both:

hetplot(f1_ad_phenotype,HPH_ad)

hetplot(f1_ad_phenotype,pHPH_ad)

We map as dominance:

hetplot(f1_d_phenotype,HPH_ad)

hetplot(f1_d_phenotype,pHPH_ad)