library(bio3d)
library(ncdf)


pdb1 <- read.pdb("1XCK_chainA_noWAT_noH.pdb")
pdb2 <- read.pdb("1SX4_A.pdb_flsq.pdb")
pdb3 <- read.pdb("1SX4_N.pdb_flsq.pdb")
pdb4 <- read.pdb("1SX3_A.pdb_flsq.pdb")


inds1 <- atom.select(pdb1, "///7:135///CA/")
inds1b <- atom.select(pdb1, "///409:518///CA/")
inds1$atom <- c(inds1$atom, inds1b$atom)
inds1$xyz <- c(inds1$xyz, inds1b$xyz)

inds2 <- atom.select(pdb2, "///1:135///CA/")
inds2b <- atom.select(pdb2, "///409:525///CA/")
inds2$atom <- c(inds2$atom, inds2b$atom)
inds2$xyz <- c(inds2$xyz, inds2b$xyz)

inds3 <- atom.select(pdb3, "///1:135///CA/")
inds3b <- atom.select(pdb3, "///409:525///CA/")
inds3$atom <- c(inds3$atom, inds3b$atom)
inds3$xyz <- c(inds3$xyz, inds3b$xyz)

inds4 <- atom.select(pdb4, "///1:135///CA/")
inds4b <- atom.select(pdb4, "///409:525///CA/")
inds4$atom <- c(inds4$atom, inds4b$atom)
inds4$xyz <- c(inds4$xyz, inds4b$xyz)


rmsd(pdb1$xyz, pdb2$xyz, a.inds=inds1$xyz, b.inds=inds2$xyz, fit=TRUE)
rmsd(pdb1$xyz, pdb3$xyz, a.inds=inds1$xyz, b.inds=inds3$xyz, fit=TRUE)
rmsd(pdb1$xyz, pdb4$xyz, a.inds=inds1$xyz, b.inds=inds4$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

t <- seq(1,5000,20)


# 128: 50ns_noWAT_1000frames_noH_chain_1.nc
# 129: 50ns_noWAT_5000frames_noH_chain_
## monomers
for ( i in 1:7 ) {
  tmptrj = read.ncdf(paste(prefix116, "traj_monomer/50ns_noWAT_5000frames_noH_chain_", i, ".nc", sep=""))
  sim116$trj = rbind(sim116$trj, tmptrj[t,])
}

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

for ( i in 1:7 ) {
  tmptrj = read.ncdf(paste(prefix171, "traj_monomer/50ns_noWAT_5000frames_noH_chain_", i, ".nc", sep=""))
  sim171$trj = rbind(sim171$trj, tmptrj[t,])
}

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


t2 <- seq(1, length(t)*7, 1)

sim116$rmsd.apo <- rmsd(pdb1$xyz, sim116$trj[t2,], a.inds=inds1$xyz, b.inds=inds1$xyz, fit=TRUE)
sim116$rmsd.trans <- rmsd(pdb3$xyz, sim116$trj[t2,], a.inds=inds3$xyz, b.inds=inds1$xyz, fit=TRUE)
sim116$rmsd.cis <- rmsd(pdb4$xyz, sim116$trj[t2,], a.inds=inds4$xyz, b.inds=inds1$xyz, fit=TRUE)
sim116$rmsd.cis2 <- rmsd(pdb2$xyz, sim116$trj[t2,], a.inds=inds2$xyz, b.inds=inds1$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$xyz, b.inds=inds$xyz, fit=TRUE)
sim129$rmsd.r <- rmsd(pdb3$xyz, sim129$trj[t2,], a.inds=ca.inds.2c7e$xyz, b.inds=inds$xyz, fit=TRUE)

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$xyz, b.inds=inds$xyz, fit=TRUE)
sim129b$rmsd.r <- rmsd(pdb3$xyz, sim129b$trj, a.inds=ca.inds.2c7e$xyz, b.inds=inds$xyz, fit=TRUE)


sim128$rmsd.r <- rmsd(pdb3$xyz, sim128$trj, a.inds=ca.inds.2c7e$xyz, b.inds=inds$xyz, fit=TRUE)

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=ca.inds.2c7e$xyz, b.inds=inds$xyz, fit=TRUE)

pdf("rmsd_R_histogram.pdf", w=15, h=5)
par(mfrow=c(1,3))
hist(sim116$rmsd.r, xlim=c(3.5, 6.5), ylim=c(0,200))
hist(sim129$rmsd.r, xlim=c(3.5, 6.5), ylim=c(0,200))
hist(sim171$rmsd.r, xlim=c(3.5, 6.5), ylim=c(0,200))
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=12)
par(mfcol=c(3,1))

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=7, byrow=T)
plot(colMeans(tmp), type='l', col="blue", ylim=c(3.5, 6))
for ( i in 1:nrow(tmp) ) {
	lines(tmp[i,], col="grey50", lty=2)
}

tmp <- matrix(sim129$rmsd.r, nrow=7, 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=7, byrow=T)
lines(colMeans(tmp), col="black")

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

dev.off()