library(bio3d)
library(gclus)
library(maptools)


load("xrayPCA.RData")


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


rd <- rmsd(xyz[, gaps.pos$f.inds])
dis <- as.dist(rd)
hc <- hclust(dis)
hc1 <- reorder.hclust(hc, dis)
id <- substr(basename(pdbs$id), 1, 6)
grps <- cutree(hc1, k=3)


prefix <- "50ns_noWAT_1000frames_noH_chain_"
trj.inds <- seq(1, 1000, by=2)


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

for ( i in 1:7 ) {
  tmptrj <- read.ncdf(paste(prefix116, prefix, i , ".nc", sep=""))
  sim116$trj = rbind(sim116$trj, tmptrj[trj.inds,])
}

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

for ( i in 1:7 ) {
  tmptrj <- read.ncdf(paste(prefix129, prefix, i , ".nc", sep=""))
  sim129$trj$atp = rbind(sim129$trj$atp, tmptrj[trj.inds,])
}

sim171 <- NULL
prefix171 <- "/net/lutefisk/slars/groel_md/1XCK_rhodanese/171_1XCK_rhodanese_ATP/results/traj_monomer/"
prefix <- "100ns_noWAT_2000frames_noATP_noH_chain_"


for ( i in 1:7 ) {
  tmptrj <- read.ncdf(paste(prefix171, prefix, i , ".nc", sep=""))
  sim171$trj$atp = rbind(sim171$trj$atp, tmptrj[trj.inds,])
}


sim117 <- NULL
prefix117 <- "/net/gulrotkake/slars/groel_md/1XCK/117_1XCK_hepta/results/traj_monomer/"
prefix <- "50ns_noWAT_1000frames_noH_chain_"

for ( i in 1:7 ) {
  tmptrj <- read.ncdf(paste(prefix117, prefix, i , ".nc", sep=""))
  sim117$trj = rbind(sim117$trj, tmptrj[trj.inds,])
}

sim128 <- NULL
prefix128 <- "/net/gulrotkake/slars/groel_md/1XCK/128_1XCK_hepta_MGATP/results/traj_monomer/"
prefix <- "50ns_noWAT_1000frames_noH_chain_"

for ( i in 1:7 ) {
  tmptrj <- read.ncdf(paste(prefix128, prefix, i , ".nc", sep=""))
  sim128$trj = rbind(sim128$trj, tmptrj[trj.inds,])
}


xyz_traj1 <- fit.xyz(fixed = pdbs$xyz[1, gaps.pos$f.inds], mobile =  sim116$trj[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, full.pdbs = FALSE, het2atom = TRUE)
xyz_traj2 <- fit.xyz(fixed = pdbs$xyz[1, gaps.pos$f.inds], mobile =  sim129$trj$atp[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, full.pdbs = FALSE, het2atom = TRUE)
xyz_traj3 <- fit.xyz(fixed = pdbs$xyz[1, gaps.pos$f.inds], mobile =  sim171$trj$atp[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, full.pdbs = FALSE, het2atom = TRUE)
xyz_traj4 <- fit.xyz(fixed = pdbs$xyz[1, gaps.pos$f.inds], mobile =  sim117$trj[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, full.pdbs = FALSE, het2atom = TRUE)
xyz_traj5 <- fit.xyz(fixed = pdbs$xyz[1, gaps.pos$f.inds], mobile =  sim128$trj[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, full.pdbs = FALSE, het2atom = TRUE)

d1 <- pca.project(xyz_traj1, pc.xray)
d2 <- pca.project(xyz_traj2, pc.xray)
d3 <- pca.project(xyz_traj3, pc.xray)
d4 <- pca.project(xyz_traj4, pc.xray)
d5 <- pca.project(xyz_traj5, pc.xray)

save(d1,d2,d3, file="projections.RData")


highlight.ind2 <-  c(grep("1XCK_A", id), grep("1SVT_A", id), grep("2C7E_A", id))


hl.col=c("orange", "green", "red")
cex.main=1.2



xlim=c(-150,300); ylim=c(-50,200);
pdf("confplot.pdf", w=12, h=3)
par(mfrow=c(1,5))
conf.plot( d1, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=xlim, ylim=ylim, ylab="PC2", xlab="PC1", mtext="A", mtext2="", main="Closed unbound", cex.main=1 )
conf.plot( d2, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=xlim, ylim=ylim, ylab="PC2", xlab="PC1", mtext="B", mtext2="", main="Closed ATP bound", cex.main=1 )
conf.plot( d3, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=xlim, ylim=ylim, ylab="PC2", xlab="PC1", mtext="C", mtext2="", main="Closed ATP bound with rhodanese", cex.main=1 )
conf.plot( d4, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=xlim, ylim=ylim, ylab="PC2", xlab="PC1", mtext="D", mtext2="", main="Closed ATP bound with rhodanese", cex.main=1 )
conf.plot( d5, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=xlim, ylim=ylim, ylab="PC2", xlab="PC1", mtext="E", mtext2="", main="Closed ATP bound with rhodanese", cex.main=1 )

dev.off()





conf.plot <- function( projection, pc.xray, projection.secondary=NULL,
                      pc.inds = c(1,2),
                      xlim=NULL, ylim=NULL, main = "", xlab="PC1", ylab="PC2",
                      hl.col=NULL, hl.ind=NULL, ablines=NULL,
                      mtext=NULL, mtext2=NULL, ... ) {

  
  if( is.null(projection.secondary) ) {
    projection.secondary = projection
    col.sec = densCols(projection[,pc.inds[1]], projection[,pc.inds[2]])
    projection = NULL
  }
  else {
    col.sec="grey90"
  }
    
  plot(projection.secondary[,pc.inds[1]], projection.secondary[,pc.inds[2]], type="p",
       main = main, xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim,
       col=col.sec,
       pch=20, ...)

  if( !is.null(projection) ) {
    points( projection[,pc.inds[1]], projection[,pc.inds[2]],
           col=densCols(projection[,pc.inds[1]], projection[,pc.inds[2]]),
           pch=20)
  }
    
  points( pc.xray$z[,pc.inds[1]], pc.xray$z[,pc.inds[2]], col = "black", pch = 16, cex = 1 )
  points( pc.xray$z[,pc.inds[1]], pc.xray$z[,pc.inds[2]], col="gray", pch = 16, cex = 0.4 )
  points( pc.xray$z[hl.ind, pc.inds[1]], pc.xray$z[hl.ind, pc.inds[2]], col = hl.col, pch = 16 )
  
  mtext(mtext, side=3, at=xlim[1]-50, cex=0.9, line=0.6)
  #mtext(mtext2, side=3, at=xlim[1]*1.2, cex=0.6, line=-1)
  text(xlim[1], ylim[2]*0.95, mtext2, pos=4)

  #abline(v=mean(projection.secondary[,pc.inds[1]]), col="grey80", lty=2)

}