library(maptools)
library(rparallel)
library(bio3d)
library(ncdf)
source("dmat_funs.R")

##pdb <- read.pdb("1XCK_reference.pdb")
pdb <- read.pdb("chainB_noWAT_noH.pdb")
pdb.ref <- pdb


seq <- seq.pdb(pdb)
s <- array(seq)
sse <- dssp(pdb)
sse <- stride(pdb)


ca.inds <- atom.select(pdb, "calpha")
dim <- length(ca.inds$atom)

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


sim180 <- NULL
prefix180 <- "/net/gulrotkake/slars/groel_md/1SX4/180_1SX4/results/traj_monomer/"

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

for ( i in 8:14 ) {
  tmptrj <- read.ncdf(paste(prefix180, prefix, i , ".nc", sep=""))
  sim180$trj$trans = rbind(sim180$trj$trans, tmptrj[trj.inds,])
}


sim181 <- NULL
prefix181 <- "/net/gulrotkake/slars/groel_md/1SX4/181_1SX4_atp/results/traj_monomer/"


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

for ( i in 8:14 ) {
  tmptrj <- read.ncdf(paste(prefix181, prefix, i , ".nc", sep=""))
  sim181$trj$trans = rbind(sim181$trj$trans, tmptrj[trj.inds,])
}


sim180$dm$cis <- dm.xyz.trj3( sim180$trj$cis, pdb,  threads=2 )
sim181$dm$cis <- dm.xyz.trj3( sim181$trj$cis, pdb, threads=2 )

sim180$dm$trans <- dm.xyz.trj3( sim180$trj$trans, pdb,  threads=2 )
sim181$dm$trans <- dm.xyz.trj3( sim181$trj$trans, pdb, threads=2 )
#save(sim180, sim181, file="sim180_181.RData")
load("sim180_181.RData")


dmat <- NULL
dmat$cis <- dcm(sim180$dm$cis, sim181$dm$cis, occupancy=0.5, diff.cut=0.5)
dmat$trans <- dcm(sim180$dm$trans, sim181$dm$trans, occupancy=0.5, diff.cut=0.5)

md.avg.diff <- NULL
md.avg.diff$cis <- sim180$dm$cis$dmat.avg - sim181$dm$cis$dmat.avg
md.avg.diff$trans <- sim180$dm$trans$dmat.avg - sim181$dm$trans$dmat.avg 

md.avg.diff=md.avg.diff$cis
dmat$current=dmat$cis
md.avg.diff=md.avg.diff$trans
dmat$current=dmat$trans

inds <- which(dmat$current != 0, arr.ind=T)

diff <- which((md.avg.diff!=0) & (dmat$current!=0), arr.ind=T )
#diff <- which((md.avg.diff!=0), arr.ind=T )
diff = cbind(diff, round(md.avg.diff[diff],2))
diff = data.frame(diff)
diff = cbind(diff,  paste( s[diff[,"row"]], ((diff[,"row"]-1)%%524)+2, sep="" ) )
diff = cbind(diff,  paste( s[diff[,"col"]], ((diff[,"col"]-1)%%524)+2, sep="" ) )
colnames(diff)<-c("row", "col", "avg.diff", "res1", "res2")
diff.table1 <- diff[abs(order(diff$avg.diff)),]

diff.table1[,c("res1", "res2", "avg.diff")]


head(diff.table1, n=20)
tail(diff.table1, n=20)


a <- paste(s[inds[,1]], ((inds[,1]-1)%%524)+2, sep="")
b <- paste(s[inds[,2]], ((inds[,2]-1)%%524)+2, sep="")
#a <- paste(s[inds[,1]], (inds[,1])+1, sep="")
#b <- paste(s[inds[,2]], (inds[,2])+1, sep="")
labels <- paste(a, b, sep="-")


pdf("180_181_TRANS_dcm.pdf")
##pdf("180_181_CIS_dcm.pdf")
plot.dcm(dmat$current, pdb, xlim=c(1,524), ylim=c(1,524),
         sse=sse, sse.grid=TRUE,
         )
#pointLabel(inds[,1], inds[,2], labels, cex=0.3, offset=10)


rows=c(1:5)
tmp1=tail(diff.table1, n=20)
pointLabel(tmp1[,"row"], tmp1[,"col"],
           paste(tmp1[,"res1"],tmp1[,"res2"],sep="-"),
           cex=0.5, offset=10)
tmp2=head(diff.table1, n=20)
pointLabel(tmp2[,"row"], tmp2[,"col"],
           paste(tmp2[,"res1"],tmp2[,"res2"],sep="-"),
           cex=0.5, offset=10)

abline(v=c(133,190,377,409), lty=2, col="grey50")
abline(h=c(133,190,377,409), lty=2, col="grey50")

dev.off()


## HOLO contacts
tmp1=tail(diff.table1, n=20)
c1 <- paste(tail(tmp1[,"row"],n=16), collapse="+")
c2 <- paste(tail(tmp1[,"col"],n=16), collapse="+")
paste(c1,c2,sep="+")

## APO contacts
tmp2=head(diff.table1, n=20)
c1 <- paste(head(tmp2[,"row"],n=13), collapse="+")
c2 <- paste(head(tmp2[,"col"],n=13), collapse="+")
paste(c1,c2,sep="+")