library(bio3d)
load("xrayPCA.RData")
source("confplot_funs.R")
library(car)

#load("xrayPCA_Tstates.RData")


pdb.closed <- read.pdb("1XCK_chainA_noWAT_noH.pdb")
ca.inds <- atom.select(pdb.closed, "calpha")

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


sim116 <- NULL
prefix116 <- "/net/gulrotkake/slars/groel_md/1XCK/116_1XCK_apo/results/traj_monomer/"
prefix <- "20-50ns_noWAT_1500frames_noH_chain_"

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

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

sim129 <- NULL
prefix129 <- "/net/gulrotkake/slars/groel_md/1XCK/129_1XCK_MGATP/results/traj_monomer/"
prefix <- "20-50ns_noWAT_1500frames_noATP_noH_chain_"

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

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


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


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


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


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

sim116$xyzfit <- NULL
sim129$xyzfit <- NULL

save(sim116, sim129, file="projected.RData")







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

conf.plot( sim129$cis$proj, pc.xray, mtext="ATP-bound (cis)", 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=-55)

conf.plot( sim129$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("B", side=3, line=1, at=-55)

conf.plot( sim116$cis$proj, pc.xray, mtext="APO (cis)", 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("C", side=3, line=1, at=-55)

conf.plot( sim116$trans$proj, pc.xray, mtext="APO (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("D", side=3, line=1, at=-55)

dev.off()









xlim=c(-60,80)
ylim=c(-50,40)
pdf("mdplot_all.pdf", width=10)
par(mfcol=c(2,2))

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

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



conf.plot( d1, pc.xray, mtext="APO", xlim=xlim, ylim=ylim,
          hl.ind=hl.ind, hl.col=hl.col)
abline(v=mean(d1[,1]), lty=2, col="grey60")
abline(h=mean(d1[,2]), lty=2, col="grey60")
mtext("C", side=3, line=1, at=-20)


conf.plot( d2, pc.xray, mtext="APO", xlim=xlim, ylim=ylim,
          hl.ind=hl.ind, hl.col=hl.col)
abline(v=mean(d2[,1]), lty=2, col="grey60")
abline(h=mean(d2[,2]), lty=2, col="grey60")
mtext("D", side=3, line=1, at=-20)


dev.off()


}