Post date: Oct 14, 2019 10:11:9 PM
I played some with visualizing the hypercube from the all forms of epistasis model. The code is in plotFitnessLandscape.R within /uufs/chpc.utah.edu/common/home/u6000989/projects/timema_genarch/color_exp_2019_gbs/Gemma_mel
library(igraph)
# 3 binary
G<-cbind(rep(c(0,1),each=4),rep(c(0,0,1,1),2),rep(c(0,1),4))
nds<-paste("s",1:8,sep="")
colnames(nds)<-"id"
jns<-matrix(NA,nrow=28,ncol=3)
jns<-as.data.frame(jns)
n<-1
for(i in 1:7){for(j in (i+1):8){
jns[n,1]<-nds[i]
jns[n,2]<-nds[j]
jns[n,3]<-sum(abs(G[i,]-G[j,]))
n<-n+1
}}
colnames(jns)<-c("from","to","weight")
jns2<-jns[which(jns[,3]==1),]
net<-graph_from_data_frame(d=jns2,vertices=nds,directed=FALSE)
l<-layout_with_gem(net)
plot(net,layout=l)
# 5 binary
2^5
#[1] 32
G<-cbind(rep(c(0,1),each=16),
rep(c(rep(0,8),rep(1,8)),2),
rep(c(rep(0,4),rep(1,4)),4),
rep(c(0,0,1,1),8),
rep(c(0,1),16))
nds<-paste("s",1:32,sep="")
nds<-as.data.frame(nds)
colnames(nds)<-"id"
(32*31)/2
jns<-matrix(NA,nrow=496,ncol=3)
jns<-as.data.frame(jns)
n<-1
for(i in 1:31){for(j in (i+1):32){
jns[n,1]<-as.character(nds[i,1])
jns[n,2]<-as.character(nds[j,1])
jns[n,3]<-sum(abs(G[i,]-G[j,]))
n<-n+1
}}
colnames(jns)<-c("from","to","weight")
jns2<-jns[which(jns[,3]==1),]
net<-graph_from_data_frame(d=jns2,vertices=nds,directed=FALSE)
l<-layout_with_gem(net)
plot(net,layout=l)
# 5, 3 genotype
3^5
#[1] 243
G<-cbind(rep(c(0,1,2),each=81),
rep(c(rep(0,27),rep(1,27),rep(2,27)),3),
rep(c(rep(0,9),rep(1,9),rep(2,9)),9),
rep(c(0,0,0,1,1,1,2,2,2),27),
rep(c(0,1,2),81))
nds<-rep(NA,243)
for(i in 1:243){
nds[i]<-paste(G[i,],collapse="")
}
nds<-as.data.frame(nds)
colnames(nds)<-"id"
(243*242)/2
jns<-matrix(NA,nrow=29403,ncol=3)
jns<-as.data.frame(jns)
n<-1
for(i in 1:242){for(j in (i+1):243){
jns[n,1]<-as.character(nds[i,1])
jns[n,2]<-as.character(nds[j,1])
jns[n,3]<-sum(abs(G[i,]-G[j,]))
n<-n+1
}}
colnames(jns)<-c("from","to","weight")
jns2<-jns[which(jns[,3]==1),]
net<-graph_from_data_frame(d=jns2,vertices=nds,directed=FALSE)
l<-layout_with_gem(net)
plot(net,layout=l)
plot(net,layout=l,vertex.size=2,vertex.color=c(rep("orange",100),rep("yellow",143)))
efGac<-scan("bmavAC.txt")
cGac<-matrix(c(-1.74489571,-1.87624877,-1.19373817,-1.27585320,-0.20343575,-0.05078484,
-0.7448957,-0.8762488,-0.1937382,-0.2758532,0.7965642,0.9492152,
0.2551043,0.1237512,0.8062618,0.7241468,1.7965642,1.9492152),nrow=3,byrow=TRUE)
exfit<-rep(NA,243)
for(i in 1:243){
gg<-rep(NA,5)
for(j in 1:5){
gg[j]<-cGac[G[i,j]+1,j]
}
y<-1
form<-y~.^5
mod<-model.matrix(form,dat=as.data.frame(t(gg)))
exfit[i]<-sum(mod[-1] * efGac)
}
exfitAC<-exfit
efGmm<-scan("bmavMM.txt")
cGmm<-matrix(c(-1.74272546,-1.83116023,-1.20107060,-1.38443852,-0.15445727,-0.04479259,
-0.7427255,-0.8311602,-0.2010706,-0.3844385,0.8455427,0.9552074,
0.2572745,0.1688398,0.7989294,0.6155615,1.8455427,1.9552074),nrow=3,byrow=TRUE)
exfit<-rep(NA,243)
for(i in 1:243){
gg<-rep(NA,5)
for(j in 1:5){
gg[j]<-cGmm[G[i,j]+1,j]
}
y<-1
form<-y~.^5
mod<-model.matrix(form,dat=as.data.frame(t(gg)))
exfit[i]<-sum(mod[-1] * efGmm)
}
exfitMM<-exfit
library(RColorBrewer)
pdf("chumashFitLands.pdf",width=22,height=11)
cs<-rev(brewer.pal(n=9,"RdYlGn"))
par(mfrow=c(1,2))
csfit<-rep(cs[9],243)
brk<-c(1,.5,0.25,0.1,-0.1,-.25,-0.5,-1)
for(i in 1:8){
a<-which(exfitAC<brk[i])
csfit[a]<-cs[i]
}
plot(net,layout=l,vertex.size=3,vertex.color=csfit)
title(main="(A) A/C fitness landscape",cex.main=1.7)
csfit<-rep(cs[9],243)
brk<-c(1,.5,0.25,0.1,-0.1,-.25,-0.5,-1)
for(i in 1:8){
a<-which(exfitMM<brk[i])
csfit[a]<-cs[i]
}
plot(net,layout=l,vertex.size=3,vertex.color=csfit)
title(main="(B) MM fitness landscape",cex.main=1.7)
dev.off()