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

##t <- seq(1,5000,20)

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, 5)

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, 5)

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, 5)
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)

pdf("rmsd_R_histogram.pdf", w=10, h=10)
par(mfrow=c(3,2))
hist(sim116$rmsd.r, xlim=c(6., 7.8))
hist(sim116$rmsd.r.2, xlim=c(6., 7.8))

hist(sim129$rmsd.r, xlim=c(6., 7.8))
hist(sim171$rmsd.r, xlim=c(6., 7.8))
hist(sim129b$rmsd.r, xlim=c(6., 7.8))
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()