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

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


load("pca_dimers.RData")
id <- substr(basename(pdbs$id), 1, 7)

col <- seq(1, length(id))
col = 0
t.inds <- c(grep("1XCK", id), grep("2NWC", id))
atp.trans <- c(grep("1SX4", id), grep("1AON", id), grep("1SVT", id), grep("1PF9", id))
atp.cis <- c(grep("1SX3", id), grep("1KP8", id))
r.inds <- grep("2C7E", id)

col.inds <- seq(1, length(id))*0
col.inds[t.inds] = 1
col.inds[atp.trans] = 2
col.inds[atp.cis] = 3
col.inds[r.inds] = 4

col[t.inds] = "red"
col[atp.trans] = "orange"
col[atp.cis] = "green"
col[r.inds] = "blue"

                    
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

d <- matrix(sim116$cis$proj[,1], ncol=1500, byrow=TRUE)
e <- matrix(sim129$cis$proj[,1], ncol=1500, byrow=TRUE)

pdf("projectionsVsTime.pdf", w=12, h=18)
par(mfrow=c(3,1))
plot(d[1,], type='l')
lines(e[1,], lty=2, col='red')
plot(d[2,], type='l')
lines(e[1,], lty=2, col='red')
plot(d[3,], type='l')
lines(e[1,], lty=2, col='red')
dev.off()


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

hl.ind <-  c(grep("1XCK_AB", id), grep("2C7E_AB", id), 
             grep("1SX3_AB", id), grep("1AON_MN", id) )
hl.col <- c("red", "blue", "green", "orange")


xlim=c(-60,255)
ylim=c(-50,90)


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_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("(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)


plot(pc.xray$z[,1], pc.xray$z[,2], col=col, 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)

radius=c(2,2,2,2,2)
if(FALSE) {
## T conf
for ( i in 1:3 ) {
  r2.inds <- which(col.inds==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(pc.xray$z[t.inds[1],1], pc.xray$z[t.inds[1],2]+10, "GroEL (apo)", cex=.9, pos=4)
text(pc.xray$z[atp.cis[1],1]+20, pc.xray$z[atp.cis[1],2], "GroEL-ATP (cis)", cex=.9, pos=4)
text(pc.xray$z[atp.trans[1],1], pc.xray$z[atp.trans[1],2], "GroEL-ADP (trans)", cex=.9, pos=4)
text(pc.xray$z[r.inds[1],1]-45, pc.xray$z[r.inds[1],2]+7, "R", cex=.9, pos=4)
#GroEL-ATP (R-state)

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.RData")
if(FALSE) {
conf.plot( sim116b$cis$proj, pc.xray, mtext="", xlim=xlim, ylim=ylim,
            hl.ind=hl.ind, hl.col=hl.col)
  abline(v=mean(sim116b$cis$proj[,1]), lty=2, col="grey60")
  abline(h=mean(sim116b$cis$proj[,2]), lty=2, col="grey60")
label="(A) GroEL (cis)"
mtext(label, side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)

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

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

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

load("projected_1XCK-rhod.RData")
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("(C) 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("(D) 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("(E) 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("(F) GroEL"[WT], " "[E434K], " ATP"[7], " (trans)"))
mtext(label, side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)


plot(pc.xray$z[,1], pc.xray$z[,2], col=col, 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)

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(pc.xray$z[t.inds[1],1], pc.xray$z[t.inds[1],2]+10, "GroEL (apo)", cex=.9, pos=4)
text(pc.xray$z[atp.cis[1],1]+20, pc.xray$z[atp.cis[1],2], "GroEL-ATP (cis)", cex=.9, pos=4)
text(pc.xray$z[atp.trans[1],1], pc.xray$z[atp.trans[1],2], "GroEL-ADP (trans)", cex=.9, pos=4)
text(pc.xray$z[r.inds[1],1]-45, pc.xray$z[r.inds[1],2]+7, "R", cex=.9, pos=4)


dev.off()






##ylim=ylim, xlim=xlim, 

pdf("xray_confplot3.pdf")
plot(pc.xray$z[,1], pc.xray$z[,2], col=col, xlab = "PC1", ylab = "PC2", 
	cex=0.8, pch=20)
text(pc.xray$z[hl.ind, 1], pc.xray$z[hl.ind, 2], id[hl.ind], cex=0.6, pos=1, col="grey50")
##text(pc.xray$z[, 1], pc.xray$z[, 2], id, cex=0.8, pos=c(4,1,2,1,4,4,4,1), col="grey50")
abline(h=0, lty=2, col="grey60")
abline(v=0, lty=2, col="grey60")

radius=c(4,4,3,3,2)

## T conf
for ( i in 1:4 ) {
  r2.inds <- which(col.inds==i)
  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(pc.xray$z[t.inds[1],1], pc.xray$z[t.inds[1],2]+5, "GroEL (apo)", cex=.9, pos=1)
text(pc.xray$z[atp.cis[1],1]+20, pc.xray$z[atp.cis[1],2], "GroEL-ATP (cis)", cex=.9, pos=4)
text(pc.xray$z[atp.trans[1],1], pc.xray$z[atp.trans[1],2], "GroEL-ADP (trans)", cex=.9, pos=4)
text(pc.xray$z[r.inds[1],1]-15, pc.xray$z[r.inds[1],2]+2, "R", cex=.9, pos=4)

dev.off()