simple epistasis model

epistasis_alleles.R

clr = c("#d53e4f","#fc8d59","#99d594")
n=500e3
hsq = 0.5
seeds = 10

par(mfrow=c(1,4))
for ( maf in c(0,0.05,0.2,0.5) ) {
  plot(0,0,xaxt="n",las=1,type="n",pch=19,ylab="Phenotype",xlab="Genotype",ylim=c(-1,30),xlim=c(0,2),bty="n",main="")
  axis(side=1,at=c(0,1,2),labels = c("AA","Aa","aa"))
  hsq_add = rep(NA,seeds)
  for ( s in 1:seeds ) {
    # generate epistasis panels
    if ( maf > 0 ){
      X1 = rbinom(n,2,maf)
      X2 = rbinom(n,2,maf)
      y = scale((X1)*(X2))*sqrt(hsq) + rnorm(n,0,sqrt(1-hsq))
    # generate left-most additive 5% freq panel 
    } else {
      X1 = rbinom(n,2,0.05)
      X2 = rbinom(n,2,0.05)
      y = scale(X1 + X2)*sqrt(hsq) + rnorm(n,0,sqrt(1-hsq))
    }
    
    # estimate the fraction of epistasis that can be captured by additivity
    hsq_add[s] = summary(lm(y ~ X1 + X2))$adj.r.sq / summary(lm(y ~ X1*X2))$adj.r.sq
    lines(c(0,1,2),c(mean(y[X1==0 & X2 == 0]),mean(y[X1==1 & X2 == 0]),mean(y[X1==2 & X2 == 0])),type="o",pch=19,col=paste(clr[1],"50",sep=''))
    lines(c(0,1,2),c(mean(y[X1==0 & X2 == 1]),mean(y[X1==1 & X2 == 1]),mean(y[X1==2 & X2 == 1])),type="o",pch=19,col=paste(clr[2],"50",sep=''))
    lines(c(0,1,2),c(mean(y[X1==0 & X2 == 2]),mean(y[X1==1 & X2 == 2]),mean(y[X1==2 & X2 == 2])),type="o",pch=19,col=paste(clr[3],"50",sep=''))
    
    }
  if ( maf == 0 ) {
    title("No epistasis: 5% frequency\n(100% additive)")
  } else {
    title(paste("Epistasis: ",maf*100,"% frequency\n(",round(mean(hsq_add)*100),"% additive)",sep=''))
  }
  
  cat( mean(hsq_add) , '\n' )
  legend("topleft",pch=19,col=clr,legend=c("BB","Bb","bb"),bty="n")
}