library(gclus)
library(maptools)
library(bio3d)
source("confplot_funs.R")
#load("xrayPCA.RData")
load("xrayPCA_Tstates.RData")



      

pdb.e434k <- read.pdb("E434K_apo_chainA_noWAT_noH.pdb")
ca.inds.e434k <- atom.select(pdb.e434k, "calpha")

trj.inds <- seq(1, 1500, by=1)

sim191 <- NULL
prefix191 <- "/net/lutefisk/slars/groel_md/1XCK/191_1XCK_E434K/results/traj_monomer/"
prefix <- "20-50ns_noWAT_1500frames_noATP_noH_chain_"

for ( i in 1:7 ) {
  tmptrj <- read.ncdf(paste(prefix191, prefix, i , ".nc", sep=""))
  sim191$xyzfit[[i]] <- fit.xyz(fixed = pdbs$xyz[1,gaps.pos$f.inds],
                                mobile = tmptrj[trj.inds, ca.inds.e434k$xyz],
                                fixed.inds = core$c1A.xyz,
                                mobile.inds = core$c1A.xyz,
                                full.pdbs = FALSE)
}

for ( i in 8:14 ) {
  tmptrj <- read.ncdf(paste(prefix191, prefix, i , ".nc", sep=""))
  sim191$xyzfit[[i]] <- fit.xyz(fixed = pdbs$xyz[1,gaps.pos$f.inds],
                                mobile = tmptrj[trj.inds, ca.inds.e434k$xyz],
                                fixed.inds = core$c1A.xyz,
                                mobile.inds = core$c1A.xyz,
                                full.pdbs = FALSE)
}

## E434K
tmptrj <- NULL
for ( i in 1:7 ) {
  tmptrj=rbind(tmptrj, sim191$xyzfit[[i]])
}
sim191$cis$proj <- pca.project(tmptrj, pc.xray)

tmptrj <- NULL
for ( i in 8:14 ) {
  tmptrj=rbind(tmptrj, sim191$xyzfit[[i]])
}
sim191$trans$proj <- pca.project(tmptrj, pc.xray)


sim191$xyzfit <- NULL
save(sim191, file="projected_1XCK-E434K.RData")







load("xrayPCA_Tstates.RData")
load("projected_1XCK-E434K.RData")
id <- substr(basename(pdbs$id), 1, 6)
hl.ind <-  c(grep("1XCK_A", id), 
             grep("1SX3_A", id), grep("1AON_N", id) )
hl.col <- c("red", "green", "orange")


at.offset=3
mtext.cex=0.75
line=.4
xlim=c(-60,120)
ylim=c(-50,60)
pdf("e434k_mdplot_all.pdf", w=6, h=4)
par(mfcol=c(2,3), mar=c(3,3,1.5,1), mgp=c(1.7,.8,0))
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")
mtext("(A) E434K-ATP (cis)", 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")
mtext("(B) E434K-ATP (trans)", side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)

load("../13_xrayPCA_EQ/xrayPCA_EQ.RData")
load("../13_xrayPCA_EQ/projected_sim162.RData")

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


hl.ind <-  c(grep("1XCK_A", id), grep("1AON_A", id),
             grep("1SX3_A", id), grep("1AON_N", id) )
hl.col <- c("red", "yellow", "green", "orange")
xlim=c(-30,23); ylim=c(-20,15);

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")
mtext("(C) E434K-ATP (cis)", 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")
mtext("(D) E434K-ATP (trans)", side=3, line=line, cex=mtext.cex, at=xlim[1]-at.offset, adj=0)

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

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

dev.off()








if (FALSE) {
conf.plot( sim161$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")
mtext("A", side=3, line=1, at=-20)

conf.plot( sim161$trans$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")
  mtext("A", side=3, line=1, at=-20)
}





xlim=c(-60,120)
ylim=c(-50,60)
pdf("e461k_mdplot_all.pdf", w=8, h=8)
par(mfcol=c(2,2))


conf.plot( sim151$cis$proj, pc.xray, mtext="ATP-free (trans)", xlim=xlim, ylim=ylim,
            hl.ind=hl.ind, hl.col=hl.col)
  abline(v=mean(proj[,1]), lty=2, col="grey60")
  abline(h=mean(proj[,2]), lty=2, col="grey60")
  mtext("A", side=3, line=1, at=-20)

conf.plot( sim151$trans$proj, pc.xray, mtext="ATP-free (trans)", xlim=xlim, ylim=ylim,
            hl.ind=hl.ind, hl.col=hl.col)
  abline(v=mean(proj[,1]), lty=2, col="grey60")
  abline(h=mean(proj[,2]), lty=2, col="grey60")
  mtext("A", side=3, line=1, at=-20)

conf.plot( sim150$cis$proj, pc.xray, mtext="ATP-free (trans)", xlim=xlim, ylim=ylim,
            hl.ind=hl.ind, hl.col=hl.col)
  abline(v=mean(proj[,1]), lty=2, col="grey60")
  abline(h=mean(proj[,2]), lty=2, col="grey60")
mtext("A", side=3, line=1, at=-20)

conf.plot( sim150$trans$proj, pc.xray, mtext="ATP-free (trans)", xlim=xlim, ylim=ylim,
            hl.ind=hl.ind, hl.col=hl.col)
  abline(v=mean(proj[,1]), lty=2, col="grey60")
  abline(h=mean(proj[,2]), lty=2, col="grey60")
  mtext("A", side=3, line=1, at=-20)

dev.off()