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

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

seq <- seq.pdb(pdb)
blast <- blast.pdb(seq)

##print(blast)
#pdf("blastplot.pdf", onefile=TRUE)
hits <- plot.blast(blast, cutoff=500 ) 

exclude.inds = grep("1AON_[HIJKLMN]|1SX4_[HIJKLMN]|1PF9_[HIJKLMN]|1SVT_[HIJKLMN]|1XCK|1KP8|1SX3|1GRU[HIJKLMN]|1SS8|2C7C[HIJKLMN]|2NWC[HIJKLMN]|3FBH|2C7E[ABCDEFG]", 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,gaps.pos$f.inds], mobile = pdbs, fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", outpath = "core_fitlsq_with2c7e/", full.pdbs = TRUE, het2atom = TRUE)

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

pc.xray <- pca.xyz(xyz[,gaps.pos$f.inds])

library(gclus)
library(maptools)

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)
fit<-kmeans( cbind(pc.xray$z[,1], pc.xray$z[,2]), 3, 20)


pdf("xray_PCA.pdf")
plot(pc.xray, cex=0.5)
dev.off()


highlight.ind <-  c(grep("1SX4_N", id),grep("1PF9_N", id),
                    grep("1AON_K", id),grep("1SVT_I", id),
                    grep("1SVT_K", id),grep("1AON_J", id),
                    
                    grep("1SX3_A", id),grep("1KP8_B", id),
                    grep("1KP8_F", id),grep("1SX3_I", id),
                    grep("1SX3_C", id),grep("2C7C_H", id),

                    grep("2C7E_A", id),grep("2C7E_G", id),
                    
                    grep("1SS8_C", id),grep("3FBH_N", id),
                    grep("1XCK_B", id),grep("1XCK_D", id),
                    grep("1SS8_A", id),grep("1XCK_N", id) )


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

save(fit, grps, highlight.ind, pc.xray, gaps.pos, core, aln, pdbs, pdb, hits,
     file="xrayPCA_Tstates_with2C7E.RData")
#load("xrayPCA_Tstates.RData")


pdf("xray_confplot4.pdf", w=6, h=6)
par(mgp=c(2,1,0))
plot(pc.xray$z[,1], pc.xray$z[,2], col=fit$cluster, xlab = "PC1", ylab = "PC2", 
        ylim=c(-35,23), xlim=c(-40,200), cex=0.8, pch=20, cex.axis=.8, cex.lab=.8)
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(2,2,2,2,2)

## 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(-8, -10, "GroEL (apo)", cex=.9)
text(18,12, "GroEL-ATP (cis)", cex=.9)
text(-25,7, "GroEL-ADP (trans)", cex=.9)
dev.off()




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

trj.inds <- seq(1, 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)
}

sim171 <- NULL
prefix171 <- "/net/lutefisk/slars/groel_md/1XCK_rhodanese/171_1XCK_rhodanese_ATP/results/traj_monomer/"
prefix <- "100ns_noWAT_2000frames_noATP_noH_chain_"
trj.inds <- seq(500, 2000, by=3)

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



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



xlim=c(-60,100)
ylim=c(-50,40)
pdf("mdplot_all.pdf", width=10, h=20)
par(mfcol=c(7,4))

for ( i in 1:7 ) {
  proj <- pca.project(sim129$xyzfit[[i]], pc.xray)
  
  conf.plot( proj, pc.xray, mtext="ATP (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=-20)
}

for ( i in 8:14 ) {
  proj <- pca.project(sim129$xyzfit[[i]], pc.xray)
  
  conf.plot( proj, pc.xray, mtext="ATP (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)
}


for ( i in 1:7 ) {
  proj <- pca.project(sim116$xyzfit[[i]], pc.xray)
  
  conf.plot( proj, pc.xray, mtext="ATP-free (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=-20)
}


for ( i in 8:14 ) {
  proj <- pca.project(sim116$xyzfit[[i]], pc.xray)
  
  conf.plot( 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()




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)

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


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


sim116$xyzfit <- NULL
sim129$xyzfit <- NULL
sim171$xyzfit <- NULL
save(sim116, sim129, file="projected.RData")







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