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

library(gridBase)
library(plotrix)
library(car)
library(cluster)

#load("xrayPCA_EQ.RData")
#source("confplot_funs.R")

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



seq <- seq.pdb(pdb)
blast <- blast.pdb(seq)
hits <- plot.blast(blast, cutoff=500 ) 
unq.ids <- unique(substr(hits$pdb.id, 1, 4))

for ( i in grep("1WF4", hits$pdb.id) ) {
  hmm <- hits$pdb.id[i]
  hits$pdb.id[i] = substr(hmm, 1, 6)
  }
                                        
files <- paste("raw_pdbs_new/split_chain/", hits$pdb.id, ".pdb", sep = "")


aln <- pdbaln(files)
pdbs <- read.fasta.pdb(aln, "", "", het2atom = T)
gaps.res <- gap.inspect(pdbs$ali)
exclude.inds <- which(gaps.res$row > 9)
##save(exclude.inds, gaps.res, seq, blast, hits, aln, pdbs, files, file="inital_aln.RData")
load("inital_aln.RData")


exclude.inds = c(grep("1WF4", hits$pdb.id),
  grep("1WE3", hits$pdb.id), grep("3C9V", hits$pdb.id),
  grep("1IOK", hits$pdb.id), grep("2CGT", hits$pdb.id),
  grep("2CGT", hits$pdb.id), grep("1GR5", hits$pdb.id),
  grep("1GRL", hits$pdb.id), grep("3CAU", hits$pdb.id) )


files.excl <- files[-exclude.inds]
aln <- pdbaln(files.excl)
pdbs <- read.fasta.pdb(aln, "", "", het2atom = T)


core <- core.find(pdbs)
xyz <- fit.xyz(fixed = pdbs$xyz[1, ], mobile = pdbs, fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", outpath = "core_fitlsq_new/", full.pdbs = TRUE, het2atom = TRUE)



gaps.pos <- gap.inspect(pdbs$xyz)
gaps.res <- gap.inspect(pdbs$ali)



pdb.closed <- read.pdb("1XCK_chainA_noWAT_noH.pdb")
##inds <- atom.select(pdb, "///500:515///CA/")
ca.inds <- atom.select(pdb.closed, "calpha")


e1 <- gaps.pos$f.inds[atom2xyz(1:134)] ## from Ala3
e2 <- gaps.pos$f.inds[atom2xyz(408:522)] ## to Leu524
eq <- c(e1,e2)

e1 <- atom2xyz(2:135)
e2 <- atom2xyz(409:523)
eq.md <- c(e1,e2)


pc.xray <- pca.xyz(xyz[,eq])
#plot(pc.xray, pch = 1 )

rd <- rmsd(xyz[,eq])
dis <- as.dist(rd)
hc <- hclust(dis)
hc1 <- reorder.hclust(hc, dis)
id <- substr(basename(pdbs$id), 1, 6)
grps <- cutree(hc1, k=4)
fit<-kmeans( cbind(pc.xray$z[,1], pc.xray$z[,2]), 4, 10)


highlight.ind <-  c(
                    grep("1GRU_M", id), grep("1SX4_N", id),
                    grep("1SVT_H", id),
                         
                    grep("1KP8_B", id),grep("1SS8_B", id),
                    grep("1XCK_B", id),grep("3FBH_A", id),
                    grep("2NWC_I", id),

                    grep("1SX4_A", id), grep("1SVT_D", id), 
                    grep("1AON_A", id),

                    grep("1SX3_B", id), grep("1SX3_N", id),
                    grep("3E76_H", id), grep("1KP8_A", id),

                    grep("2C7E_A", id),

                    grep("2C7C_K", id),
                    grep("2C7C_J", id), grep("2C7C_H", id),
                    grep("2C7C_L", id), grep("2C7C_M", id),
                    grep("2C7C_N", id), grep("2C7C_B", id), 
                    grep("2C7C_A", id), grep("2C7C_F", id),
                    
                    grep("2C7D_I", id),
                    grep("2C7D_H", id), grep("2C7D_A", id),
                    grep("2C7D_L", id), grep("2C7D_M", id) )


pos=c(
  4,4,4,
  4,4,2,4,4,
  2,1,3,
  4,4,4,4,
  1,
  4,4,4,4,4,4,4,4,4,
  4,4,4,4
  )

save(pc.xray, highlight.ind, pos, hits, gaps.pos, grps, fit, files, aln, pdbs, blast, core, xyz, file="xrayPCA_EQ.RData")
load(file="xrayPCA_EQ.RData")

plot.pca(pc.xray)


load("colorsForClustering.RData")


pdf("xray_confplot3.pdf")
plot(pc.xray$z[,1], pc.xray$z[,2], col=mycolors, xlab = "PC1", ylab = "PC2", 
	ylim=c(-10,13), xlim=c(-27,28), cex=0.8, pch=20)
text(pc.xray$z[highlight.ind, 1], pc.xray$z[highlight.ind, 2], id[highlight.ind], cex=0.6, pos=pos, 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)

cluster <- seq(1, length(id)) * 0
cluster[which(mycolors=="red")] <- 1
cluster[which(mycolors=="green")] <- 2
cluster[which(mycolors=="orange")] <- 3
cluster[which(mycolors=="yellow")] <- 4

## T conf
for ( i in 1:4 ) {
  r2.inds <- which(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(1, 9, "GroEL-ATP (cis)", cex=.9)
text(-6.5,-8, "GroEL (apo)", cex=.9)
text(-20,7, "GroEL-ADP (trans)", cex=.9)
text(18,1.5, "GroEL-ADP (cis)", cex=.9)
dev.off()












prefix <- "50ns_noWAT_5000frames_noH_chain_"
trj.inds <- seq(2500, 5000, by=5)


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$cis = rbind(sim116$trj$cis, tmptrj[trj.inds,])
}
for ( i in 8:14 ) {
  tmptrj <- read.ncdf(paste(prefix116, prefix, i , ".nc", sep=""))
  sim116$trj$trans = rbind(sim116$trj$trans, tmptrj[trj.inds,])
}

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

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

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


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

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

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

if (FALSE) {
sim180 <- NULL
prefix180 <- "/net/gulrotkake/slars/groel_md/1SX4/180_1SX4/results/traj_monomer/"

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

for ( i in 8:14 ) {
  tmptrj <- read.ncdf(paste(prefix180, prefix, i , ".nc", sep=""))
  sim180$trj$free = rbind(sim180$trj$free, tmptrj[trj.inds,])
}

sim181 <- NULL
prefix181 <- "/net/gulrotkake/slars/groel_md/1SX4/181_1SX4_atp/results/traj_monomer/"

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

for ( i in 8:14 ) {
  tmptrj <- read.ncdf(paste(prefix181, prefix, i , ".nc", sep=""))
  sim181$trj$atp = rbind(sim181$trj$atp, tmptrj[trj.inds,])
}
}







xyz_traj1 <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  sim116$trj$cis[, 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, ], mobile =  sim116$trj$trans[, 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, ], mobile =  sim129$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, ], mobile =  sim129$trj$free[, 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, ], mobile =  sim171$trj$atp[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, full.pdbs = FALSE, het2atom = TRUE)
xyz_traj6 <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  sim171$trj$free[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, full.pdbs = FALSE, het2atom = TRUE)

if (FALSE) {
xyz.180 <- NULL
xyz.181 <- NULL
xyz.180$adp <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  sim180$trj$adp[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, full.pdbs = FALSE, het2atom = TRUE)
xyz.180$free <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  sim180$trj$free[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, full.pdbs = FALSE, het2atom = TRUE)

xyz.181$adp <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  sim181$trj$adp[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, full.pdbs = FALSE, het2atom = TRUE)
xyz.181$atp <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  sim181$trj$atp[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, full.pdbs = FALSE, het2atom = TRUE)
}

d1 <- pca.project(xyz_traj1[,eq.md], pc.xray)
d2 <- pca.project(xyz_traj2[,eq.md], pc.xray)
d3 <- pca.project(xyz_traj3[,eq.md], pc.xray)
d4 <- pca.project(xyz_traj4[,eq.md], pc.xray)
d5 <- pca.project(xyz_traj5[,eq.md], pc.xray)
d6 <- pca.project(xyz_traj6[,eq.md], pc.xray)

if (FALSE) {
d180a <- pca.project(xyz.180$adp[,eq.md], pc.xray)
d180b <- pca.project(xyz.180$free[,eq.md], pc.xray)
d181a <- pca.project(xyz.181$adp[,eq.md], pc.xray)
d181b <- pca.project(xyz.181$atp[,eq.md], pc.xray)
}


save(d1,d2,d3,d4,d5,d6, file="projected.RData")




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

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


sim161 <- NULL
prefix161 <- "/net/lutefisk/slars/groel_md/mutants/161_E434K_apo/results/traj_monomer/"
prefix <- "50ns_noWAT_1000frames_noH_chain_"

for ( i in 1:7 ) {
  tmptrj <- read.ncdf(paste(prefix161, prefix, i , ".nc", sep=""))
  sim161$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(prefix161, prefix, i , ".nc", sep=""))
  sim161$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)
}


sim162 <- NULL
prefix162 <- "/net/lutefisk/slars/groel_md/mutants/162_E434K_MGATP/results/traj_monomer/"
prefix <- "50ns_noWAT_1000frames_noH_chain_"


for ( i in 1:7 ) {
  tmptrj <- read.ncdf(paste(prefix162, prefix, i , ".nc", sep=""))
  sim162$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(prefix162, prefix, i , ".nc", sep=""))
  sim162$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, sim161$xyzfit[[i]])
}
sim161$cis$proj <- pca.project(tmptrj[,eq.md], pc.xray)

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

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

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

sim162$xyzfit <- NULL
sim161$xyzfit <- NULL
save(sim161, sim162, file="projected_sim162.RData")




hl.ind <-  c(grep("1XCK_A", id), grep("1SX4_A", id),
             grep("1SX3_A", id), grep("1SX4_N", id) )
hl.col <- c("orange", "green", "yellow", "red")

pdf("mdplot_all.pdf", width=10)
par(mfcol=c(2,3))

conf.plot( d3, pc.xray, mtext="ATP (cis)", xlim=c(-20,25), ylim=c(-20,20),
           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("A", side=3, line=1, at=-20)

conf.plot( d4, pc.xray, mtext="ATP free (trans)", xlim=c(-20,25), ylim=c(-20,20),
          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("B", side=3, line=1, at=-20)

conf.plot( d5, pc.xray, mtext="ATP (cis)", xlim=c(-20,25), ylim=c(-20,20),
          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("C", side=3, line=1, at=-20)

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

conf.plot( d1, pc.xray, mtext="APO", xlim=c(-20,25), ylim=c(-20,20),
          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("E", side=3, line=1, at=-20)

conf.plot( d2, pc.xray, mtext="APO", xlim=c(-20,25), ylim=c(-20,20),
          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("E", side=3, line=1, at=-20)

dev.off()


pdf("mdplot_1sx4.pdf", width=8)
par(mfcol=c(2,2))

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

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

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

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

dev.off()






#plot(pc.xray$z[,2], pc.xray$z[,3], col=grps+1, xlab = "PC2", ylab = "PC3", cex=0.8, pch=20)
#text(pc.xray$z[highlight.ind, 2], pc.xray$z[highlight.ind, 3], id[highlight.ind], 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")



#plot(pc.xray$z[,1], pc.xray$z[,3], col=grps+1, xlab = "PC2", ylab = "PC3", cex=0.8, pch=20)
#text(pc.xray$z[highlight.ind, 1], pc.xray$z[highlight.ind, 3], id[highlight.ind], 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")



vi <- mktrj.pca(pc.xray, pc=1, file="pc1_eq.pdb")
vi <- mktrj.pca(pc.xray, pc=2, file="pc2_eq.pdb", mag=4, step=0.50)
vi <- mktrj.pca(pc.xray, pc=3, file="pc3_eq.pdb", mag=4, step=0.50)
vi <- mktrj.pca(pc.xray, pc=4, file="pc4_eq.pdb", mag=4, step=0.50)
vi <- mktrj.pca(pc.xray, pc=5, file="pc5_eq.pdb", mag=4, step=0.50)
#vi <- mktrj.pca(pc.xray, pc=2, file="pc2_eq.pdb")