library(bio3d)
load("xrayPCA.RData")

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


source("overlap_funs.R")
k <- grep("1XCK_A", pdbs$id)
l <- grep("1SVT_A", pdbs$id)
dv <- diff.vector(xyz, c(k,l),  gaps.pos$f.inds)



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


sim116 <- NULL
prefix116 <- "/net/gulrotkake/slars/groel_md/1XCK/116_1XCK_apo/results/traj_monomer/"

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

trj.inds <- seq(1, 5000, by=50)
sim121 <- NULL
prefix121 <- "/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/"
prefix <- "50ns_noWAT_noATP_5000frames_noH_p"

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


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

sim121$fit.xyz <- fit.xyz(fixed = pdbs$xyz[1, gaps.pos$f.inds], mobile =  sim121$trj[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, full.pdbs = FALSE, het2atom = TRUE)


sim116$pc <- pca.xyz(sim116$fit.xyz)
sim129$atp$pc <- pca.xyz(sim129$atp$fit.xyz)
#sim129$free$pc <- pca.xyz(sim129$free$fit.xyz)
sim121$pc <- pca.xyz(sim121$fit.xyz)


for ( i in 1:7 ) {
  end=i*length(trj.inds)
  start=end-length(trj.inds)+1
  sim116$A$pc$indv.units[[i]] <- pca.xyz(sim116$fit.xyz[start:end,])
}
#for ( i in 8:14 ) {
#  end=i*length(trj.inds)
#  start=end-length(trj.inds)+1
#  sim116$B$pc$indv.units[[i]] <- pca.xyz(sim116$fit.xyz[start:end,])
#}

for ( i in 1:7 ) {
  end=i*length(trj.inds)
  start=end-length(trj.inds)+1
  pcs <- pca.xyz(sim129$atp$fit.xyz[start:end,])
  sim129$atp$pc$indv.units[[i]] = pcs
}

#for ( i in 1:7 ) {
#  end=i*length(trj.inds)
#  start=end-length(trj.inds)+1
#  pcs <- pca.xyz(sim171$atp$fit.xyz[start:end,])
#  sim171$atp$pc$indv.units[[i]] = pcs
#}




#a <- c()
#for ( i in 1:7 ) {
#  r <- rmsip(pc.xray, sim129$atp$pc$indv.units[[i]], subset=10)
#  a = c(a, r$rmsip)
#}

#b <- c()
#for ( i in 1:7 ) {
#  r <- rmsip(pc.xray, sim116$A$pc$indv.units[[i]], subset=10)
#  b = c(b, r$rmsip)
#}


#pdf("rmsip_plot.pdf", w=10, h=4)
#par(mfcol=c(2,4), mar=c(1.,1,0.5,0.5))
#for ( i in 1:7 ) {
#  r <- rmsip(pc.xray, sim116$A$pc$indv.units[[i]], subset=10)
#  image(1:10, 1:10, r$matrix, col=gray(10:0/10), zlim=c(0,1),
#        xlab="MD-PCs", ylab="X-ray PCs")
#             #zlim=c(0,1), 
#}
#dev.off()




cex=1.2
pdf("RMSIP_and_Overlap.pdf", w=9, h=2)
layout(matrix(c(1,2,3,4), 1, 4, byrow=T),
       widths=c(2,2,2,3), heights=c(2,2,2,2))
par(mgp=c(1.5,0.5,0), mar=c(3,3,1,1))

r <- rmsip(sim116$pc, sim121$pc, subset=10)
image(1:10, 1:10, r$matrix, col=gray(20:0/20), zlim=c(0,1),
      xlab="ATP-free MD PCs", ylab="Monomer PCs")
mtext(3, at=-1.5, line=-.6, text="A", cex=cex)
mtext(paste("RMSIP:", round(r$rmsip,2)), side=3, line=-1.5, at=1, cex=0.75, adj=0)

r <- rmsip(pc.xray, sim129$atp$pc, subset=10)
image(1:10, 1:10, r$matrix, col=gray(20:0/20), zlim=c(0,1),
      xlab="ATP-bound MD PCs", ylab="X-ray PCs")
mtext(3, at=-1.5, line=-.6, text="B", cex=cex)
mtext(paste("RMSIP:", round(r$rmsip,2)), side=3, line=-1.5, at=1, cex=0.75, adj=0)

r <- rmsip(sim116$pc, sim129$atp$pc, subset=10)
image(1:10, 1:10, r$matrix, col=gray(20:0/20), zlim=c(0,1),
      xlab="ATP-bound MD PCs", ylab="ATP-free MD PCs")
mtext(3, at=-1.5, line=-.6, text="C", cex=cex)
mtext(paste("RMSIP:", round(r$rmsip,2)), side=3, line=-1.5, at=1, cex=0.75, adj=0)

a <- c()
b <- c()
c <- c()
d <- c()
plot(-1, ylim=c(0.1,0.9), xlim=c(0,20), ylab="Cumulative overlap", xlab="Mode index")
mtext(3, at=-2.5, line=-.6, text="D", cex=cex)
for ( i in 1:7 ) {
  o <- overlap(sim116$A$pc$indv.units[[i]], dv)
  lines(o$c, lty=2, col="grey50")
  a=rbind(a,o$c)
}
#for ( i in 8:14 ) {
#  o <- overlap(sim116$B$pc$indv.units[[i]], dv)
#  lines(o$c, lty=2, col="grey50")
#  d=rbind(d,o$c)
#}

for ( i in 1:7 ) {
  o <- overlap(sim129$atp$pc$indv.units[[i]], dv)
  lines(o$c, lty=3, col="grey50")
  b=rbind(b,o$c)
}

#for ( i in 1:7 ) {
#  o <- overlap(sim171$atp$pc$indv.units[[i]], dv)
#  lines(o$c, lty=4, col="grey50")
#  c=rbind(c,o$c)
#}

lines(colMeans(a), lty=1, type="o", cex=0.5,col="blue")
#lines(colMeans(d), lty=1, type="o", cex=0.5,col="blue")
lines(colMeans(b), lty=1, type="o", cex=0.5,col="red")
#lines(colMeans(c), lty=1, type="o", cex=0.5,col="orange")

sig <- c()
for ( i in 1:20 ) {
  t <- t.test(a[,i], b[,i])
  print(t$p.value)
  if ( t$p.value < 0.05 ){
    sig=c(sig, i)
    #print(t$p.value)
  }
}
points(sig, colMeans(a)[sig], col="darkgreen", pch=16)

dev.off()




plot.pca(sim116$pc)
plot.pca(sim129$atp$pc)
plot.pca(sim129$free$pc)
dev.off()
a <- mktrj.pca(sim116$pc, pc=1, file="116_pc1.pdb", mag=1.5, step=0.175)
a <- mktrj.pca(sim116$pc, pc=2, file="116_pc2.pdb", mag=1.5, step=0.175)
a <- mktrj.pca(sim116$pc, pc=3, file="116_pc3.pdb", mag=1.5, step=0.175)
a <- mktrj.pca(sim116$pc, pc=4, file="116_pc4.pdb", mag=1.5, step=0.175)
a <- mktrj.pca(sim116$pc, pc=5, file="116_pc5.pdb", mag=1.5, step=0.175)
a <- mktrj.pca(sim116$pc, pc=7, file="116_pc7.pdb", mag=1.5, step=0.175)

a <- mktrj.pca(sim129$atp$pc, pc=1, file="129_atp_pc1.pdb", mag=1.5, step=0.175)
a <- mktrj.pca(sim129$atp$pc, pc=2, file="129_atp_pc2.pdb", mag=1.5, step=0.175)
a <- mktrj.pca(sim129$atp$pc, pc=3, file="129_atp_pc3.pdb", mag=1.5, step=0.175)
a <- mktrj.pca(sim129$atp$pc, pc=6, file="129_atp_pc6.pdb", mag=1.5, step=0.175)







source("overlap_funs.R")
k <- grep("1XCK_A", pdbs$id)
l <- grep("1SVT_A", pdbs$id)
dv <- diff.vector(xyz, c(k,l),  gaps.pos$f.inds)
sim116$pc$overlap <- overlap(sim116$pc, dv)
sim116$pc$rmsip <- rmsip(pc.xray, sim116$pc)

sim129$atp$pc$overlap <- overlap(sim129$atp$pc, dv)
sim129$atp$pc$rmsip <- rmsip(pc.xray, sim129$atp$pc)


pdf("overlaps.pdf")
par(mfcol=c(2,1))
plot(sim116$pc$overlap$v)
lines(sim129$atp$pc$overlap$v, lty=2, col=2)

plot(sim116$pc$overlap$c)
lines(sim129$atp$pc$overlap$c, lty=2, col=2)
dev.off()




plot(o$v, type='o', ylim=c(0,0.7))
lines(o2$v, lty=2, col="gray50")

plot(o$c, type='o')
lines(o2$c, lty=2, col="gray50")
dev.off()