library(bio3d)
library(car)
source("confplot_funs.R")

pdb <- read.pdb("1XCK_A.pdb")

load("xrayPCA.RData")
id <- substr(basename(pdbs$id), 1, 6)


highlight.ind <-  c(grep("1SX4_N", id),grep("1PF9_N", id),
                    grep("1AON_K", id),grep("1SVT_I", id),
                    grep("1SVT_K", id),grep("1AON_M", id),
                    
                    grep("1SX3_A", id),grep("1KP8_B", id),
                    grep("1KP8_F", id),grep("1SX3_G", id),
                    grep("1SX3_C", id),grep("2C7C_H", id),
                    
                    grep("1SS8_C", id),grep("3FBH_N", id),
                    grep("1XCK_B", id),grep("1XCK_D", id),
                    grep("1SS8_A", id),grep("1XCK_N", id) )


apo <- grep("1XCK|3FBH|2NWC_[HIJKLMN]|1SS8", id)
cis <- grep("1KP8|1SX3", id)
trans <- grep("1AON_[HIJKLMN]|1SX4_[HIJKLMN]|1PF9_[HIJKLMN]|1SVT_[HIJKLMN]", id)
r <- grep("2C7E_[ABCDEFG]", id)

colors = rep(0, length(id))
colors[apo] = "red"
colors[cis] = "green"
colors[trans] = "orange"
colors[r] = "blue"
                    
pos=c(
  4,4,4,4,4,1,
  4,4,4,1,4,4,
  4,4,4,4,4,4
  )

cex.axis=1.1
cex.lab=1.1
cex=0.9
at.offset=3
mtext.cex=0.68
line=.4


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

load("projected.RData")

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

hl.ind <-  c(grep("1XCK_A", id), grep("2C7E_A", id), 
             grep("1SX3_A", id), grep("1AON_N", id) )
hl.col <- c("red", "blue", "green", "orange")


xlim=c(-60,200)
ylim=c(-35,100)


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_1XCK-rhod.RData")


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

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




load("projected_mutant_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("(G) 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("(H) GroEL"[E434K], " 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("(I) 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("(J) GroEL"[WT], " "[E434K], " ATP"[7], " (trans)"))
mtext(label, side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)


load("xrayPCA.RData")

xlim=c(-60,200)
ylim=c(-35,40)

plot(pc.xray$z[,1], pc.xray$z[,2], col=colors, 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("(K) X-ray", side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)

radius=c(2,2,2,2,2)

if(FALSE) {
## T conf
for ( i in 1:3 ) {
  r2.inds <- which(fit$cluster==i)
  #r2.inds=r2.inds[-r.inds]
  v1 <- as.vector(pc.xray$z[r2.inds,1])
  v2 <- as.vector(pc.xray$z[r2.inds,2])
  m <- matrix(c(v1,v2), ncol=2)
  ellipse(c(mean(v1),mean(v2)), cov(m), radius[i], lty="dashed", col="grey80", center.pch=0, lwd=0.7)
  #text(mean(v1)-5, mean(v2)+5, "Cis/Trans t", cex=0.9, col="black")
}
}

text(-10, 0, "GroEL (apo)", cex=.9, pos=4)
text(10,25, "GroEL-ATP (cis)", cex=.9, pos=4)
text(-50,-30, "GroEL-ADP (trans)", cex=.9, pos=4)


dev.off()