library(bio3d)
#library(ncdf)

pdb1 <- read.pdb("1XCK_heptamer_noH.pdb")
pdb2 <- read.pdb("1SX4_cis_noWAT_MMPBSA.noH.pdb.1")
pdb3 <- read.pdb("2C7E_heptamer.pdb")

inds <- atom.select(pdb1, "///1:3668///CA/")
inds.2 <- atom.select(pdb2, "///1:3668///CA/")
inds.3 <- atom.select(pdb3, "///2:525///CA/")

rmsd(pdb1$xyz, pdb3$xyz, a.inds=inds$xyz, b.inds=inds.3$xyz, fit=TRUE)
rmsd(pdb1$xyz, pdb2$xyz, a.inds=inds$xyz, b.inds=inds.2$xyz, fit=TRUE)


prefix116 <- "/net/gulrotkake/slars/groel_md/1XCK/116_1XCK_apo/results/"
prefix129 <- "/net/gulrotkake/slars/groel_md/1XCK/129_1XCK_MGATP/results/"
prefix117 <- "/net/gulrotkake/slars/groel_md/1XCK/117_1XCK_hepta/results/"
prefix128 <- "/net/gulrotkake/slars/groel_md/1XCK/128_1XCK_hepta_MGATP/results/"
prefix171 <- "/net/lutefisk/slars/groel_md/1XCK_rhodanese/171_1XCK_rhodanese_ATP/results/"
prefix129b <- "/net/lutefisk/slars/groel_md/1XCK/129_1XCK_MGATP_temp350/results/"

sim116 <- NULL
sim129 <- NULL
sim129b <- NULL
sim117 <- NULL
sim128 <- NULL
sim171 <- NULL

## fetch projected data
#load("../31_pca_dimers/projected_50nsForTrajPlots.RData")
#e <- matrix(sim129$cis$proj[,1], nrow=6, byrow=TRUE)
#f <- matrix(sim171$cis$proj[,1], nrow=6, byrow=TRUE)


## read heptamer trajectories
#sim116$trj = read.ncdf(paste(prefix116, "traj_heptamer/50ns_noWAT_5000frames_noH_hept_1.nc", sep=""))
#sim116$trj2 = read.ncdf(paste(prefix116, "traj_heptamer/50ns_noWAT_5000frames_noH_hept_2.nc", sep=""))
#sim129$trj = read.ncdf(paste(prefix129, "traj_heptamer/50ns_noWAT_1000frames_noATP_noH_hept_1.nc", sep=""))
#sim129b$trj = read.ncdf(paste(prefix129b, "traj_heptamer/50ns_noWAT_1000frames_wATP_hept_1.nc",  sep=""))
#sim171$trj = read.ncdf(paste(prefix171, "traj_heptamer/50ns_noWAT_1000frames_noH_hept_1.nc", sep=""))


t2 <- seq(1, 1000, 1)
#sim116$rmsd.closed <- rmsd(pdb1$xyz, sim116$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
#sim116$rmsd.open <- rmsd(pdb2$xyz, sim116$trj[t2,], a.inds=inds.2$xyz, b.inds=inds$xyz, fit=TRUE)
#sim116$rmsd.r <- rmsd(pdb3$xyz, sim116$trj[t2,], a.inds=inds.3$xyz, b.inds=inds$xyz, fit=TRUE)

#sim116$rmsd.closed.2 <- rmsd(pdb1$xyz, sim116$trj2[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
#sim116$rmsd.r.2 <- rmsd(pdb3$xyz, sim116$trj2[t2,], a.inds=inds.3$xyz, b.inds=inds$xyz, fit=TRUE)

#sim129$rmsd.closed <- rmsd(pdb1$xyz, sim129$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
#sim129$rmsd.open <- rmsd(pdb2$xyz, sim129$trj[t2,], a.inds=inds.2$xyz, b.inds=inds$xyz, fit=TRUE)
#sim129$rmsd.r <- rmsd(pdb3$xyz, sim129$trj[t2,], a.inds=inds.3$xyz, b.inds=inds$xyz, fit=TRUE)

t2 <- seq(1, 800, 1)
#sim129b$rmsd.closed <- rmsd(pdb1$xyz, sim129b$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
#sim129b$rmsd.open <- rmsd(pdb2$xyz, sim129b$trj[t2,], a.inds=inds.2$xyz, b.inds=inds$xyz, fit=TRUE)
#sim129b$rmsd.r <- rmsd(pdb3$xyz, sim129b$trj[t2,], a.inds=inds.3$xyz, b.inds=inds$xyz, fit=TRUE)

t2 <- seq(1, 1000, 1)
#sim171$rmsd.closed <- rmsd(pdb1$xyz, sim171$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
#sim171$rmsd.open <- rmsd(pdb2$xyz, sim171$trj[t2,], a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
#sim171$rmsd.r <- rmsd(pdb3$xyz, sim171$trj[t2,], a.inds=inds.3$xyz, b.inds=inds$xyz, fit=TRUE)

## read monomer trajectories 
## monomers
if (FALSE) {
t <- seq(1, 5000, 5)
for ( i in 1:0 ) {
  tmptrj = read.ncdf(paste(prefix116, "traj_monomer/50ns_noWAT_5000frames_noH_chain_", i, ".nc", sep=""))
  sim116$monomertrj = rbind(sim116$monomertrj, tmptrj[t,])
}

for ( i in 1:0 ) {
  tmptrj = read.ncdf(paste(prefix129, "traj_monomer/50ns_noWAT_5000frames_noH_chain_", i, ".nc", sep=""))
  sim129$monomertrj = rbind(sim129$monomertrj, tmptrj[t,])
}

t <- seq(1,4000,5)
for ( i in 1:0 ) {
  tmptrj = read.ncdf(paste(prefix129b, "traj_monomer/500ns_noWAT_5000frames_noATP_noH_chain_", i, ".nc", sep=""))
  sim129b$monomertrj = rbind(sim129b$monomertrj, tmptrj[t,])
}

t <- seq(1, 5000, 5)
for ( i in 1:0 ) {
  tmptrj = read.ncdf(paste(prefix171, "traj_monomer/50ns_noWAT_5000frames_noH_chain_", i, ".nc", sep=""))
  sim171$monomertrj = rbind(sim171$monomertrj, tmptrj[t,])
}
}
pdb.1xckA <- read.pdb('1XCK_chainA_noWAT_noH.pdb')
inds.1xckA <- atom.select(pdb.1xckA, "///9:515///CA/")
pdb.2c7eA <- read.pdb("2C7E_A.pdb_flsq.pdb")
inds.2c7eA <- atom.select(pdb.2c7eA, "///10:516///CA/")

#sim116$rmsd.r.monomer <- rmsd(pdb.2c7eA$xyz, sim116$monomertrj, a.inds=inds.2c7eA$xyz, b.inds=inds.1xckA$xyz, fit=TRUE)
#sim129$rmsd.r.monomer <- rmsd(pdb.2c7eA$xyz, sim129$monomertrj, a.inds=inds.2c7eA$xyz, b.inds=inds.1xckA$xyz, fit=TRUE)
#sim129b$rmsd.r.monomer <- rmsd(pdb.2c7eA$xyz, sim129b$monomertrj, a.inds=inds.2c7eA$xyz, b.inds=inds.1xckA$xyz, fit=TRUE)
#sim171$rmsd.r.monomer <- rmsd(pdb.2c7eA$xyz, sim171$monomertrj, a.inds=inds.2c7eA$xyz, b.inds=inds.1xckA$xyz, fit=TRUE)

#sim129$rmsd.matrix <- matrix(sim129$rmsd.r.monomer, nrow=7, byrow=T)

sim116$trj <-NULL
sim129$trj <-NULL
sim129b$trj <-NULL
sim116$trj2 <- NULL
sim171$trj <-NULL
sim116$monomertrj <- NULL
sim129$monomertrj <- NULL
sim171$monomertrj <- NULL
sim129b$monomertrj <- NULL
#save(sim116, sim129, sim129b, sim171, file="rmsd.RData")
load("rmsd.RData")

sim116$rmsd.r.monomer <- matrix(sim116$rmsd.r.monomer, nrow=7, byrow=T)
sim129$rmsd.r.monomer <- matrix(sim129$rmsd.r.monomer, nrow=7, byrow=T)

pdf('rmsd_aposimToR_means.pdf', w=12, h=4)
plot(colMeans(sim116$rmsd.r.monomer), type='l', axes=FALSE, xlab=xlab, ylab=ylab, cex.lab=cex )
lines(colMeans(sim129$rmsd.r.monomer), col='red')
dev.off()

xlab="Time (ns)"
ylab="RMSD (Å)"
at=c(0,200,400,600,800,1000)
line=0.7
at.offset=10
mtext.cex=1.1
cex=1.2


pdf('rmsd_aposimToR.pdf', w=12, h=20)
par(mfrow=c(7,1))
plot(sim116$rmsd.r.monomer[1,], type='l', axes=FALSE, xlab=xlab, ylab=ylab, cex.lab=cex )
lines(sim129$rmsd.r.monomer[1,], col='red')
abline(h=5.3, col="grey70", lty=2)
plot(sim116$rmsd.r.monomer[2,], type='l', axes=FALSE, xlab=xlab, ylab=ylab, cex.lab=cex )
lines(sim129$rmsd.r.monomer[2,], col='red')
abline(h=5.3, col="grey70", lty=2)
plot(sim116$rmsd.r.monomer[3,], type='l', axes=FALSE, xlab=xlab, ylab=ylab, cex.lab=cex )
lines(sim129$rmsd.r.monomer[3,], col='red')
abline(h=5.3, col="grey70", lty=2)
plot(sim116$rmsd.r.monomer[4,], type='l', axes=FALSE, xlab=xlab, ylab=ylab, cex.lab=cex )
lines(sim129$rmsd.r.monomer[4, ], col='red')
abline(h=5.3, col="grey70", lty=2)
plot(sim116$rmsd.r.monomer[5, ], type='l', axes=FALSE, xlab=xlab, ylab=ylab, cex.lab=cex )
lines(sim129$rmsd.r.monomer[5, ], col='red')
abline(h=5.3, col="grey70", lty=2)
plot(sim116$rmsd.r.monomer[6, ], type='l', axes=FALSE, xlab=xlab, ylab=ylab, cex.lab=cex )
lines(sim129$rmsd.r.monomer[6,], col='red')
abline(h=5.3, col="grey70", lty=2)
plot(sim116$rmsd.r.monomer[7,], type='l', axes=FALSE, xlab=xlab, ylab=ylab, cex.lab=cex )
lines(sim129$rmsd.r.monomer[7,], col='red')
abline(h=5.3, col="grey70", lty=2)
dev.off()


pdf('rmsd_all.pdf', w=12, h=14)
par(mfrow=c(4,1))

## Plotting heptameric stuff
plot(sim116$rmsd.r, type='l', col="blue", ylim=c(6, 8), , axes=FALSE, xlab=xlab, ylab=ylab, cex.lab=cex  )
axis(2, cex.axis=cex ); axis(1, at=at, labels=c(0,10,20,30,40,50), cex.axis=cex ); box();
lines(sim116$rmsd.r.2, col="blue", lty=2)
lines(sim129$rmsd.r, col="red")
#lines(sim129b$rmsd.r, col="black")
lines(sim171$rmsd.r, col="green")
abline(h=7.33, col="grey70", lty=2)

label2=expression(paste("GroEL", " ATP"[7], " (cis)"))
label3=expression(paste("GroEL"[Rhod], " ATP"[7], " (cis)"))
legend("topright", legend=c("GroEL (ATP-free)", label2, label3), col=c("blue", "red", "green"), lty=1, ncol=3, cex=cex )
mtext("(A) Cis heptameric ring", side=3, line=line, cex=mtext.cex, at=0-at.offset, adj=0)

## Plotting subunit stuff
plot(sim129$rmsd.r.monomer[4,], type='l', axes=FALSE, xlab=xlab, ylab=ylab, cex.lab=cex, col='red' )
lines(sim116$rmsd.r.monomer[5,], col='blue')
axis(2, cex.axis=cex ); axis(1, at=at, labels=c(0,10,20,30,40,50), cex.axis=cex ); box();
abline(h=5.3, col="grey70", lty=2)
#label=expression(paste("(B) Cis subunit "[7], " cis subunit"))
label="(B) Cis subunit "
mtext(label, side=3, line=line, cex=mtext.cex, at=0-at.offset, adj=0)

## Plotting heptamer at 350K 
plot(sim129b$rmsd.r, type='l', col="red", ylim=c(6, 8), , axes=FALSE, xlab=xlab, ylab=ylab, cex.lab=cex  )
axis(2, cex.axis=cex ); axis(1, at=at[1:5], labels=c(50,60,70,80,90), cex.axis=cex ); box();
abline(h=7.33, col="grey70", lty=2)
label="(C) Cis heptameric ring (350K)"
mtext(label, side=3, line=line, cex=mtext.cex, at=0-at.offset, adj=0)

## Plotting subunit at 350K
plot(sim129b$rmsd.r.monomer[2401:3200], type='l', axes=FALSE, xlab=xlab, ylab=ylab, cex.lab=cex, col='red' )
axis(2, cex.axis=cex ); axis(1, at=at[1:5], labels=c(50,60,70,80,90), cex.axis=cex ); box();
abline(h=5.3, col="grey70", lty=2)
label="(D) Cis subunit (350K)"
mtext(label, side=3, line=line, cex=mtext.cex, at=0-at.offset, adj=0)

dev.off()
df("rmsd_R_histogram.pdf", w=10, h=8)
par(mfrow=c(2,2))
hist(sim116$rmsd.r, xlim=c(6.3, 7.8), xlab=ylab, main="GroEL (ATP-free) cis ring")
hist(sim116$rmsd.r.2, xlim=c(6.3, 7.8), xlab=ylab, main="GroEL (ATP-free) trans ring")

hist(sim129$rmsd.r, xlim=c(6.3, 7.8), xlab=ylab, main="GroEL (ATP) cis ring")
hist(sim171$rmsd.r, xlim=c(6.3, 7.8), xlab=ylab, main="GroEL (Rhodanese ATP) cis ring")
#hist(sim129b$rmsd.r, xlim=c(6.3, 7.8), xlab=ylab, main="GroEL (ATP-free) cis ring"
dev.off()

w <- wilcox.test(sim116$rmsd.r, sim129$rmsd.r, paired = FALSE, alternative="greater")
t <- t.test(sim116$rmsd.r, sim129$rmsd.r, paired = FALSE, alternative="greater")

t.test(sim116$rmsd.r, sim129$rmsd.r)


sim116$trj=NULL
sim129$trj=NULL
sim171$trj=NULL

#save(sim116,sim129,sim171, file="rmsd.RData")
load("rmsd.RData")

#rmsd(pdb1$xyz, pdb2$xyz, a.inds=inds$xyz, b.inds=inds$xyz, fit=TRUE)
#rmsd(pdb1$xyz, pdb3$xyz, a.inds=inds$xyz, b.inds=ca.inds.2c7e$xyz, fit=TRUE)

max(sim116$rmsd.closed)
max(sim129$rmsd.closed)
max(sim171$rmsd.closed)

min(sim116$rmsd.r)
min(sim129$rmsd.r)
min(sim171$rmsd.r)

min(sim116$rmsd.open)
min(sim129$rmsd.open)
min(sim171$rmsd.open)



mean(sim116$rmsd.closed)
mean(sim129$rmsd.closed)
mean(sim171$rmsd.closed)

var(sim116$rmsd.r)
var(sim129$rmsd.r)
var(sim171$rmsd.r)

mean(sim116$rmsd.r)
mean(sim129$rmsd.r)
mean(sim171$rmsd.r)

mean(sim116$rmsd.open)
mean(sim129$rmsd.open)
mean(sim171$rmsd.open)


tmp <- matrix(sim116$rmsd.closed, nrow=7, byrow=T)
mean(tmp[,100:500])
var(as.vector(tmp[,100:500]))

tmp <- matrix(sim129$rmsd.closed, nrow=7, byrow=T)
mean(tmp[,100:500])
var(as.vector(tmp[,100:500]))


tmp <- matrix(sim171$rmsd.closed, nrow=7, byrow=T)
mean(tmp[,100:500])
var(as.vector(tmp[,100:500]))




pdf("rmsd_toAll.pdf", w=12, h=8)
par(mfcol=c(2,1))


plot(sim116$rmsd.closed, type='l', col="blue", ylim=c(1, 6))
lines(sim116$rmsd.closed.2, col="blue", lty=2)
lines(sim129$rmsd.closed, col="red")
lines(sim129b$rmsd.closed, col="black")
lines(sim171$rmsd.closed, col="green")


plot(sim116$rmsd.r, type='l', col="blue", ylim=c(6, 8))
lines(sim116$rmsd.r.2, col="blue", lty=2)
lines(sim129$rmsd.r, col="red")
lines(sim129b$rmsd.r, col="black")
lines(sim171$rmsd.r, col="green")


dev.off()




tmp <- matrix(sim116$rmsd.closed, nrow=7, byrow=T)
plot(colMeans(tmp), type='l', col="blue", ylim=c(1, 3.2))
tmp <- matrix(sim129$rmsd.closed, nrow=7, byrow=T)
lines(colMeans(tmp), col="red")
tmp <- matrix(sim171$rmsd.closed, nrow=7, byrow=T)
lines(colMeans(tmp), col="black")


tmp <- matrix(sim116$rmsd.r, nrow=7, byrow=T)
plot(colMeans(tmp), type='l', col="blue", ylim=c(4.5, 6))
tmp <- matrix(sim129$rmsd.r, nrow=7, byrow=T)
lines(colMeans(tmp), col="red")
tmp <- matrix(sim171$rmsd.r, nrow=7, byrow=T)
lines(colMeans(tmp), col="black")



tmp <- matrix(sim116$rmsd.open, nrow=7, byrow=T)
plot(colMeans(tmp), type='l', col="blue", ylim=c(10.5, 12.5))
tmp <- matrix(sim129$rmsd.open, nrow=7, byrow=T)
lines(colMeans(tmp), col="red")
tmp <- matrix(sim171$rmsd.open, nrow=7, byrow=T)
lines(colMeans(tmp), col="black")



dev.off()

pdf("rmsd_toR.pdf", w=12, h=6)
par(mfcol=c(1,1))
tmp <- matrix(sim116$rmsd.r, nrow=6, byrow=T)
plot(colMeans(tmp), type='l', col="blue", ylim=c(5, 8))
for ( i in 1:nrow(tmp) ) {
	lines(tmp[i,], col="grey50", lty=2)
}

tmp <- matrix(sim129$rmsd.r, nrow=6, byrow=T)
lines(colMeans(tmp), col="red")
for ( i in 1:nrow(tmp) ) {
        lines(tmp[i,], col="grey50", lty=3)
}

tmp <- matrix(sim171$rmsd.r, nrow=6, byrow=T)
##lines(colMeans(tmp), col="black")

tmp <- matrix(sim129b$rmsd.r, nrow=6, byrow=T)
lines(colMeans(tmp), col="green")

dev.off()