library(bio3d)
library(ncdf)

load("/net/lutefisk/slars/groel_md/bio3d_groel/xrayPCA.RData")

pdb.open <- read.pdb("/net/gulrotkake/slars/groel_nma_korgpil/1SVT/chainA/start.pdb")
pdb.closed <- read.pdb("1XCK_reference.pdb")
pdb <- pdb.closed
ca.inds <- atom.select(pdb.open, "calpha")
ca.inds.closed <- atom.select(pdb.closed, "calpha")


#diff <- read.table("/net/gulrotkake/slars/groel_nma_korgpil/1SVT/chainA/uhm_diffconf.dat", sep=" ")
#diff=as.numeric(diff)

diff.id <-  c(grep("1XCK_A", pdbs$id), grep("1SVT_A", pdbs$id))
diff <- as.numeric(xyz[diff.id[1],]-xyz[diff.id[2],])
diff = normalizedVector(diff[1:1572])



#trj0 <- read.ncdf("../traj/200ns_short_2000frames.nc")
trj0 <- read.ncdf("../traj/300ns_short_15000frames_noWAT.nc")


trj1 <- read.ncdf("/net/gulrotkake/slars/groel_md/1SVT/123_1SVT_chainA_apo_seed1/results/20ns_2000frames_noWAT.nc")
trj2 <- read.ncdf("/net/gulrotkake/slars/groel_md/1SVT/123_1SVT_chainA_apo_seed2/results/20ns_2000frames_noWAT.nc")
trj3 <- read.ncdf("/net/gulrotkake/slars/groel_md/1SVT/123_1SVT_chainA_apo_seed3/results/20ns_2000frames_noWAT.nc")
trj4 <- read.ncdf("/net/gulrotkake/slars/groel_md/1SVT/123_1SVT_chainA_apo_seed4/results/20ns_2000frames_noWAT.nc")
trj5 <- read.ncdf("/net/gulrotkake/slars/groel_md/1SVT/123_1SVT_chainA_apo_seed5/results/20ns_2000frames_noWAT.nc")


xyz0 <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  trj0[1:10000, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz300ns <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  trj0[1:15000, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)


pc.traj0 <- pca.xyz(xyz0)
mktrj.pca(pc.traj0, pc=1, file="pc1.pdb")

## project x-ray structures onto Trajectory PC
xray.proj <- pca.project(xyz[, gaps.pos$f.inds], pc.traj0)

## project trajectory onto Conf. difference vector
diff.proj <- sweep(xyz0, 2, pc.traj0$mean) %*% (diff)
diff.proj.6 <- sweep(xyz0[1:300,], 2, pc.traj0$mean) %*% (diff)
traj.proj.6 <- sweep(xyz0[1:300,], 2, pc.traj0$mean) %*% (pc.traj0$U)

## project x-ray structures onto Conf.difference vector
xray.proj.diff <- sweep(xyz[, gaps.pos$f.inds], 2, pc.traj0$mean) %*% (diff)





pdf("projection.pdf", h=4, w=12)
par(mfrow=c(1,3), mar=c(4,4,2,2))


plot( pc.traj0$z[,2], diff.proj, type="p", main = "", xlab="PC1", ylab="Diff conf", col=densCols(pc.traj0$z[,2], diff.proj), pch=20, ylim=c(-150, 320))

points(traj.proj.6[,2], diff.proj.6, col = "green", pch = 16 )
points(xray.proj[diff.id,2], xray.proj.diff[diff.id,], col = c("orange", "red"), pch = 20 )


plot( pc.traj0$z[,2], pc.traj0$z[,1], type="p", main = "", xlab="PC1", ylab="PC2", col=densCols(p0[,2], p0[,1]), pch=20)
points(traj.proj.6[,2], traj.proj.6[,1], col = "green", pch = 16 )
points(xray.proj[diff.id,2], xray.proj[diff.id,1], col = c("orange", "red"), pch = 20 )

plot( pc.traj0$z[,2], pc.traj0$z[,3], type="p", main = "", xlab="PC1", ylab="PC2", col=densCols(p0[,2], p0[,3]), pch=20)
points(traj.proj.6[,2], traj.proj.6[,3], col = "green", pch = 16 )
points(xray.proj[diff.id,2], xray.proj[diff.id,3], col = c("orange", "red"), pch = 20 )


dev.off()


















xyz1 <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  trj1[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz2 <- fit.xyz(fixed =  pdbs$xyz[1, ], mobile =  trj2[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, db.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz3 <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  trj3[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz4 <- fit.xyz(fixed =  pdbs$xyz[1, ], mobile =  trj4[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, db.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz5 <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  trj5[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)



#p0 <- pca.project(xyz0[, gaps.pos$f.inds], pc.xray)
#p1 <- pca.project(xyz1[, gaps.pos$f.inds], pc.xray)
#p2 <- pca.project(xyz2[, gaps.pos$f.inds], pc.xray)
#p3 <- pca.project(xyz3[, gaps.pos$f.inds], pc.xray)
#p4 <- pca.project(xyz4[, gaps.pos$f.inds], pc.xray)
#p5 <- pca.project(xyz5[, gaps.pos$f.inds], pc.xray)


pc.traj0 <- pca.xyz(xyz0)
mktrj.pca(pc.traj0, pc=1, file="pc1.pdb")


#pc.traj0$U <- rbind(pc.traj0$U, diff)
xray.proj <- pca.project(xyz[, gaps.pos$f.inds], pc.traj0)
p0 <- pca.project(xyz0[, gaps.pos$f.inds], pc.traj0)



#mean=NULL
#mean <- apply(xyz0[1:2000, gaps.pos$f.inds], 2, mean) ## mean structure
p0.diff <- sweep(xyz0[, gaps.pos$f.inds], 2, pc.traj0$mean) %*% (-diff)
xraj.diff.proj <- sweep(xyz[, gaps.pos$f.inds], 2, pc.traj0$mean) %*% (-diff)



#mean.short=NULL
#mean.short <- apply(xyz0[1:100, gaps.pos$f.inds], 2, mean) ## mean structure
p0.short.diff <- sweep(xyz0[1:100, gaps.pos$f.inds],2, pc.traj0$mean) %*% (-diff)
p0.short <- pca.project(xyz0[1:100, gaps.pos$f.inds], pc.traj0)




pdf("projections.pdf", h=3, w=9)
par(mfrow=c(1,3), mar=c(4,4,2,2))
plot( p0[,1], p0.diff, type="p", main = "", xlab="PC1", ylab="Diff conf", col=densCols(p0[,1], p0.diff), pch=20)
points(p0.short[,1], p0.short.diff, col="green", pch = 16, cex = 1 )
points(xraj.diff.proj[diff.id],  xray.proj[diff.id, 1], col = c("orange", "red", "red"), pch = 16 )
mtext("A", side=3, line=-1.7, at=-300, cex=1.3)

plot( pc.traj0$z[,1], pc.traj0$z[,2], type="p", main = "", xlab="PC1", ylab="PC2", col=densCols(p0[,1], p0[,2]), pch=20)
points(p0.short[,1], p0.short[,2], col="green", pch = 16, cex = 1 )
points(xray.proj[diff.id, 1], xray.proj[diff.id, 2], col = c("orange", "red", "red"), pch = 16 )

mtext("B", side=3, line=-1.7, at=-300, cex=1.3)

plot( p0[,1], p0[,2], type="p", main = "", xlab="PC1", ylab="PC3", col=densCols(p0[,1], p0[,2]), pch=20)
points(p0.short[,1], p0.short[,3], col="green", pch = 16, cex = 1 )
points(xray.proj[diff.id[1], 1], xray.proj[diff.id[2], 3], col = c("orange", "red", "red"), pch = 16 )
mtext("C", side=3, line=-1.7, at=-300, cex=1.3)

dev.off()




"normalizedVector" <-
  function(v) {
    return( v/sqrt(dotProduct(v,v)) )
  }

"dotProduct" <-
  function(a,b) {
    o <- colSums( matrix( a*b, length(a) ) )
    return(o)
}







#xyz0 <- fit.xyz(fixed = pdb$xyz[ca.inds.closed$xyz], mobile = trj0[, ca.inds$xyz], pdb.path = "", pdbext = "", outpath = "core_fitlsq/", full.pdbs = FALSE, het2atom = TRUE)
#xyz1 <- fit.xyz(fixed = pdb$xyz[ca.inds.closed$xyz], mobile = trj1[, ca.inds$xyz], pdb.path = "", pdbext = "", outpath = "core_fitlsq/", full.pdbs = FALSE, het2atom = TRUE)
#xyz2 <- fit.xyz(fixed = pdb$xyz[ca.inds.closed$xyz], mobile = trj2[, ca.inds$xyz], pdb.path = "", pdbext = "", outpath = "core_fitlsq/", full.pdbs = FALSE, het2atom = TRUE)
#xyz3 <- fit.xyz(fixed = pdb$xyz[ca.inds.closed$xyz], mobile = trj3[, ca.inds$xyz], pdb.path = "", pdbext = "", outpath = "core_fitlsq/", full.pdbs = FALSE, het2atom = TRUE)
#xyz4 <- fit.xyz(fixed = pdb$xyz[ca.inds.closed$xyz], mobile = trj4[, ca.inds$xyz], pdb.path = "", pdbext = "", outpath = "core_fitlsq/", full.pdbs = FALSE, het2atom = TRUE)
#xyz5 <- fit.xyz(fixed = pdb$xyz[ca.inds.closed$xyz], mobile = trj5[, ca.inds$xyz], pdb.path = "", pdbext = "", outpath = "core_fitlsq/", full.pdbs = FALSE, het2atom = TRUE)


mean=NULL
mean <- apply(xyz0, 2, mean) ## mean structure
diff.p0 <- sweep(xyz0,2,mean) %*% (-diff)
mean=NULL

mean <- apply(xyz1, 2, mean) ## mean structure
diff.p1 <- sweep(xyz1,2,mean) %*% -diff
mean=NULL

mean <- apply(xyz2, 2, mean) ## mean structure
diff.p2 <- sweep(xyz2,2,mean) %*% -diff
mean=NULL

mean <- apply(xyz3, 2, mean) ## mean structure
diff.p3 <- sweep(xyz3,2,mean) %*% -diff
mean=NULL

mean <- apply(xyz4, 2, mean) ## mean structure
diff.p4 <- sweep(xyz4,2,mean) %*% -diff
mean=NULL

mean <- apply(xyz5, 2, mean) ## mean structure
diff.p5 <- sweep(xyz5,2,mean) %*% -diff



pdf("tmp.pdf")
par(mfrow=c(2,3))

plot( p0[,1], p0[,2],  type="p",
     col=densCols(p0[,1], p0[,2]), 
     pch=20 )


plot( p1[,1], p1[,2],  type="p",
     col=densCols(p0[,1], p0[,2]), 
     pch=20 )


plot( p2[,1], p2[,2],  type="p",
     col=densCols(p0[,1], p0[,2]), 
     pch=20 )


plot( p3[,1], p3[,2],  type="p",
     col=densCols(p0[,1], p0[,2]), 
     pch=20 )


plot( p4[,1], p4[,2],  type="p",
     col=densCols(p0[,1], p0[,2]), 
     pch=20 )


plot( p5[,1], p5[,2],  type="p",
     col=densCols(p0[,1], p0[,2]), 
     pch=20 )

dev.off()














plot( p1[,1], diff.p1, type="p",
     col=densCols(p1[,1], diff.p1), 
     pch=20 )


plot( p2[,1], diff.p2, type="p",
     col=densCols(p2[,1], diff.p2), 
     pch=20 )

plot( p3[,1], diff.p3, type="p",
     col=densCols(p3[,1], diff.p3), 
     pch=20 )


plot( p4[,1], diff.p4, type="p",
     col=densCols(p4[,1], diff.p4), 
     pch=20 )


plot( p5[,1], diff.p5, type="p",
     col=densCols(p5[,1], diff.p5), 
     pch=20 )




dev.off()




pdf("ugh..pdf", h=12)
par(mfcol=c(6,1), mar=c(3,2,1,1))

mean=NULL
mean <- apply(xyz0, 2, mean) ## mean structure
z <- sweep(xyz0,2,mean) %*% diff
plot(z,  col="gray", type="p", pch = 16, cex = 0.4, ylab="Diff conf", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,50,100,150,200))
box()

mean <- apply(xyz1, 2, mean) ## mean structure
z <- sweep(xyz1,2,mean) %*% diff
plot(z,  col="gray", type="p", pch = 16, cex = 0.4, ylab="Diff conf", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,5,10,15,20))
box()

mean <- apply(xyz2, 2, mean) ## mean structure
z <- sweep(xyz2,2,mean) %*% diff
plot(z,  col="gray", type="p", pch = 16, cex = 0.4, ylab="Diff conf", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,5,10,15,20))
box()

mean <- apply(xyz3, 2, mean) ## mean structure
z <- sweep(xyz3,2,mean) %*% diff
plot(z,  col="gray", type="p", pch = 16, cex = 0.4, ylab="Diff conf", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,5,10,15,20))
box()

mean <- apply(xyz4, 2, mean) ## mean structure
z <- sweep(xyz4,2,mean) %*% diff
plot(z,  col="gray", type="p", pch = 16, cex = 0.4, ylab="Diff conf", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,5,10,15,20))
box()

mean <- apply(xyz5, 2, mean) ## mean structure
z <- sweep(xyz5,2,mean) %*% diff
plot(z,  col="gray", type="p", pch = 16, cex = 0.4, ylab="Diff conf", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,5,10,15,20))
box()

dev.off()




inds <- atom.select(pdb.closed, "///10:520///CA/")


pdf("proj_vs_rmsd0.pdf", h=5)
par(mfcol=c(2,1), mar=c(3,3,2,1))

mean <- apply(xyz0, 2, mean) ## mean structure
z <- sweep(xyz0,2,mean) %*% diff
plot(z,  col="gray", type="p", pch = 16, cex = 0.4, ylab="Diff conf", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,50,100,150,200))
box()

rmsd<- rmsd( pdb.closed$xyz, trj, a.inds=inds$xyz, b.inds=inds$xyz, fit = TRUE )
plot(rmsd,  col="gray", type="p", pch = 16, cex = 0.4, ylab="RMSD", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,50,100,150,200))
box()
dev.off()





pdf("proj_vs_rmsd1.pdf", h=5)
par(mfcol=c(2,1), mar=c(3,2,1,1))

mean <- apply(xyz1, 2, mean) ## mean structure
z <- sweep(xyz1,2,mean) %*% diff
plot(z,  col="gray", type="p", pch = 16, cex = 0.4, ylab="Diff conf", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,50,100,150,200))
box()

pdb.closed <- read.pdb("/net/gulrotkake/slars/groel_nma_korgpil/1SVT/chainA/1xck_reference_fitToStart.pdb")
ca.inds <- atom.select(pdb, "///10:520///CA/")
rmsd<- rmsd( pdb$xyz, xyz1, a.inds=inds$xyz, b.inds=c(31:1563), fit = TRUE )
plot(rmsd,  col="gray", type="p", pch = 16, cex = 0.4, ylab="rmsd", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,5,10,15,20))
box()
dev.off()





pdf("proj_vs_rmsd2.pdf", h=5)
par(mfcol=c(2,1), mar=c(3,2,1,1))

mean <- apply(xyz2, 2, mean) ## mean structure
z <- sweep(xyz2,2,mean) %*% diff
plot(z,  col="gray", type="p", pch = 16, cex = 0.4, ylab="Diff conf", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,50,100,150,200))
box()

pdb.closed <- read.pdb("/net/gulrotkake/slars/groel_nma_korgpil/1SVT/chainA/1xck_reference_fitToStart.pdb")
ca.inds <- atom.select(pdb, "///10:520///CA/")

rmsd<- rmsd( pdb$xyz, xyz2, a.inds=inds$xyz, b.inds=c(31:1563), fit = TRUE )
plot(rmsd,  col="gray", type="p", pch = 16, cex = 0.4, ylab="rmsd", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,5,10,15,20))
box()
dev.off()





pdf("proj_vs_rmsd3.pdf", h=5)
par(mfcol=c(2,1), mar=c(3,2,1,1))

mean <- apply(xyz3, 2, mean) ## mean structure
z <- sweep(xyz3,2,mean) %*% diff
plot(z,  col="gray", type="p", pch = 16, cex = 0.4, ylab="Diff conf", xlab="Time (ns)", axes=F)
axis(2)
axis(1, at=c(1,500,1000,1500,2000), label=c(1,50,100,150,200))
box()

pdb.closed <- read.pdb("/net/gulrotkake/slars/groel_nma_korgpil/1SVT/chainA/1xck_reference_fitToStart.pdb")
ca.inds <- atom.select(pdb, "///10:520///CA/")
rmsd<- rmsd( pdb$xyz, xyz3, a.inds=inds$xyz, b.inds=c(31:1563), fit = TRUE )
plot(rmsd,  col="gray", type="p", pch = 16, cex = 0.4, ylab="rmsd", xlab="Time (ns)", axes=F)
dev.off()