library(bio3d)
library(car)
source("confplot_funs.R")
load(file="xrayPCA_EQ.RData")



if(FALSE) {
id <- substr(basename(pdbs$id), 1, 7)

col <- seq(1, length(id))
col = 0
adp.cis <- grep("2C7E", id)
atp.trans <- c(grep("1SX4", id), grep("1AON", id), grep("1SVT", id), grep("1PF9", id))
atp.cis <- c(grep("1SX3_[ABCDEFG]", id), grep("1KP8", id))
apo.inds <- <- c(grep("1XCK", id), grep("2NWC", id), grep("3E76", id), grep("1OEL", id),
                 grep("1SS8", id), )
col[t.inds] = "red"
col[atp.trans] = "orange"
col[atp.cis] = "green"
col[r.inds] = "blue"
               }



id <- substr(basename(pdbs$id), 1, 6)

hl.ind <-  c(grep("1XCK_A", id), grep("1SX4_A", id),
             grep("1SX3_A", id), grep("1SX4_N", id) )
hl.col <- c("red", "yellow", "green", "orange")


e1 <- gaps.pos$f.inds[atom2xyz(1:134)] ## from Ala3
e2 <- gaps.pos$f.inds[atom2xyz(408:522)] ## to Leu524
eq <- c(e1,e2)


cex.axis=1.1
cex.lab=1.1
cex=0.9
at.offset=3
##ylim=c(-35,23); xlim=c(-42,50);
mtext.cex=0.68
line=.4


ylim=c(-15,15)
xlim=c(-28,28)

pdf("md_conformerplot_main.pdf", w=8, h=4)
par(mfrow=c(2,4), mgp=c(2,.9,0), mar=c(3,3,2,1))


load("projected_WT.RData")
conf.plot( sim116$cis$proj, pc.xray, mtext="", xlim=xlim, ylim=ylim,
            hl.ind=hl.ind, hl.col=hl.col)
  abline(v=mean(sim116$cis$proj[,1]), lty=2, col="grey60")
  abline(h=mean(sim116$cis$proj[,2]), lty=2, col="grey60")
mtext("(A) GroEL (cis)", side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)

conf.plot( sim116$trans$proj, pc.xray, mtext="", xlim=xlim, ylim=ylim,
            hl.ind=hl.ind, hl.col=hl.col)
  abline(v=mean(sim116$trans$proj[,1]), lty=2, col="grey60")
  abline(h=mean(sim116$trans$proj[,2]), lty=2, col="grey60")
 mtext("(B) GroEL (trans)", side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)


conf.plot( sim129$cis$proj, pc.xray, mtext="", xlim=xlim, ylim=ylim,
            hl.ind=hl.ind, hl.col=hl.col)
  abline(v=mean(sim129$cis$proj[,1]), lty=2, col="grey60")
  abline(h=mean(sim129$cis$proj[,2]), lty=2, col="grey60")
label=expression(paste("(C) GroEL ATP"[7], " (cis)"))
mtext(label, side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)

conf.plot( sim129$trans$proj, pc.xray, mtext="", xlim=xlim, ylim=ylim,
            hl.ind=hl.ind, hl.col=hl.col)
  abline(v=mean(sim129$trans$proj[,1]), lty=2, col="grey60")
  abline(h=mean(sim129$trans$proj[,2]), lty=2, col="grey60")
label=expression(paste("(D) GroEL ATP"[7], " (trans)"))
mtext(label, side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)



load("projected_E434K.RData")
conf.plot(sim162$cis$proj, pc.xray, mtext="", xlim=xlim, ylim=ylim,
          hl.ind=hl.ind, hl.col=hl.col)
abline(v=mean(sim162$cis$proj[,1]), lty=2, col="grey60")
abline(h=mean(sim162$cis$proj[,2]), lty=2, col="grey60")
label=expression(paste("(E) GroEL"[E434K], " ATP"[7], " (cis)"))
mtext(label, side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)

conf.plot( sim162$trans$proj, pc.xray, mtext="", xlim=xlim, ylim=ylim,
            hl.ind=hl.ind, hl.col=hl.col)
abline(v=mean(sim162$trans$proj[,1]), lty=2, col="grey60")
abline(h=mean(sim162$trans$proj[,2]), lty=2, col="grey60")
label=expression(paste("(F) GroEL"[E434K], " ATP"[7], " (trans)"))
mtext(label, side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)



load("colorsForClustering.RData")
ylim=c(-10,13);
xlim=c(-26,23);
plot(pc.xray$z[,1], pc.xray$z[,2], col=mycolors, xlab = "PC1", ylab = "PC2", 
     ylim=ylim, xlim=xlim, cex=cex, pch=20,
     cex.axis=cex.axis, cex.lab=cex.lab)
abline(h=0, lty=2, col="grey60")
abline(v=0, lty=2, col="grey60")
mtext("(G) X-ray", side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)


text(0, 9, "GroEL (apo)", cex=.9, pos=4)
text(-23,-9, "GroEL-ATP (cis)", cex=.9, pos=4)
text(-26,6, "GroEL-ADP (cis)", cex=.9, pos=4)
text(24,-0.5, "GroEL-ADP (trans)", cex=.9, pos=2)
dev.off()










pdf("md_conformerplot_supporting.pdf", w=8, h=4)
par(mfrow=c(2,4), mgp=c(2,.9,0), mar=c(3,3,2,1))
load("projected_GroEL_rhod.RData")
## 20-50ns ##
conf.plot(sim171$cis$proj, pc.xray, mtext="", xlim=xlim, ylim=ylim,
          hl.ind=hl.ind, hl.col=hl.col)
abline(v=mean(sim171$cis$proj[,1]), lty=2, col="grey60")
abline(h=mean(sim171$cis$proj[,2]), lty=2, col="grey60")
label=expression(paste("(A) GroEL"[Rhod], " ATP"[7], " (cis)"))
mtext(label, side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)

conf.plot( sim171$trans$proj, pc.xray, mtext="", xlim=xlim, ylim=ylim,
            hl.ind=hl.ind, hl.col=hl.col)
abline(v=mean(sim171$trans$proj[,1]), lty=2, col="grey60")
abline(h=mean(sim171$trans$proj[,2]), lty=2, col="grey60")
label=expression(paste("(B) GroEL"[Rhod], " ATP"[7], " (trans)"))
mtext(label, side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)


load("projected_1XCK-E434K.RData")

conf.plot(sim191$cis$proj, pc.xray, mtext="", xlim=xlim, ylim=ylim,
          hl.ind=hl.ind, hl.col=hl.col)
abline(v=mean(sim191$cis$proj[,1]), lty=2, col="grey60")
abline(h=mean(sim191$cis$proj[,2]), lty=2, col="grey60")
label=expression(paste("(C) GroEL"[WT], " "[E434K], " ATP"[7], " (cis)"))
mtext(label, side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)

conf.plot( sim191$trans$proj, pc.xray, mtext="", xlim=xlim, ylim=ylim,
            hl.ind=hl.ind, hl.col=hl.col)
abline(v=mean(sim191$trans$proj[,1]), lty=2, col="grey60")
abline(h=mean(sim191$trans$proj[,2]), lty=2, col="grey60")
label=expression(paste("(D) GroEL"[WT], " "[E434K], " ATP"[7], " (trans)"))
mtext(label, side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)


if(FALSE){
load(file="xrayPCA_EQ.RData")
e1 <- gaps.pos$f.inds[atom2xyz(1:134)] ## from Ala3
e2 <- gaps.pos$f.inds[atom2xyz(408:522)] ## to Leu524
eq <- c(e1,e2)

fit<-kmeans( cbind(pc.xray$z[,1], pc.xray$z[,2]), 4, 30)
id <- substr(basename(pdbs$id), 1, 6)

hl.ind <-  c(grep("1XCK_A", id), grep("1SX4_A", id),
             grep("1SX3_A", id), grep("1SX4_N", id) )

c1 <- fit$cluster[grep("1XCK_A", id)]
c2 <- fit$cluster[grep("1SX4_A", id)]
c3 <- fit$cluster[grep("1SX3_A", id)]
c4 <- fit$cluster[grep("1SX4_N", id)]

t1 <- which(fit$cluster==c1)
t2 <- which(fit$cluster==c2)
t3 <- which(fit$cluster==c3)
t4 <- which(fit$cluster==c4)

cluster.col=fit$cluster
cluster.col[t1]="red"
cluster.col[t2]="yellow"
cluster.col[t3]="green"
cluster.col[t4]="orange"
}

#mycolors=cluster.col
#save(mycolors, file="colorsForClustering.RData")
load("colorsForClustering.RData")


ylim=c(-10,13);
xlim=c(-26,23);
plot(pc.xray$z[,1], pc.xray$z[,2], col=mycolors, xlab = "PC1", ylab = "PC2", 
     ylim=ylim, xlim=xlim, cex=cex, pch=20,
     cex.axis=cex.axis, cex.lab=cex.lab)
abline(h=0, lty=2, col="grey60")
abline(v=0, lty=2, col="grey60")
mtext("(E) X-ray", side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)


text(0, 9, "GroEL (apo)", cex=.9, pos=4)
text(-23,-9, "GroEL-ATP (cis)", cex=.9, pos=4)
text(-26,6, "GroEL-ADP (cis)", cex=.9, pos=4)
text(24,-0.5, "GroEL-ADP (trans)", cex=.9, pos=2)
dev.off()