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

##pdb <- read.pdb("1XCK_2reference.pdb")
pdb <- read.pdb("1XCK_chainA_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(1, 5000, by=5)
frames <- length(trj.inds)


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

for ( i in 1:5 ) {
  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/"
#prefix129 <- "~/gulrotkake/groel_md/1XCK/129_1XCK_MGATP/results/traj_monomer/"


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




#wsize=frames * (2/100)
wsize=100

for ( i in 1:19 ) {

  apo <-  NULL
  holo <-  NULL

  for ( j in 1:7 ) {
    start=((j-1)*frames)+1+(wsize*(i-1))
    end=start+wsize-1

    print(start)
    print(end)
    
    apo = rbind(apo, sim116$trj[c(start:end),])
    holo = rbind(holo, sim129$trj$atp[c(start:end),])
  }

    
  dm.apo <- dm.xyz.trj3( apo, pdb,  threads=4 )
  print(dim(dm.apo$dmat))
  dm.holo <- dm.xyz.trj3( holo, pdb, threads=4 )

  dmat <- dcm(dm.apo, dm.holo, occupancy=0.5, diff.cut=0.2)
  save(dm.apo, dm.holo, dmat, file=paste("time3_", i, ".dat", sep=""))
  
  #png(paste("time3_", i ,".png", sep=""))
  #plot.dcm(dmat, pdb, xlim=c(1,524), ylim=c(1,524),
  #         sse=sse, sse.grid=TRUE )
  #dev.off()
}



for ( i in 1:19 ) {
  
  j=i
  while ( nchar(j) < 3 ) {
    j=paste("0",j, sep="")
  }

  load(paste("time3_", i, ".dat", sep=""))
  png(paste("time_", j, ".png", sep=""))
  plot.dcm(dmat, pdb, xlim=c(1,524), ylim=c(1,524),
           sse=sse, sse.grid=TRUE )
  dev.off()

}



## distance matrix only:
sse.grid=TRUE; sse.grid.col="gray80"; sse.grid.lty="dotted";
helix.col="gray20"; sheet.col="gray80";
xlab='Residue No'; ylab='Residue No';
cex=0.6; pch=c(2,20); col=c('green', 'blue');
labels.cex=NULL; labels.offset=5;
sse=NULL; top=FALSE; bot=TRUE; left=TRUE;

for ( i in 1:nrow(sim116$trj) ) {
   
  j=i
  while ( nchar(j) < 3 ) {
    j=paste("0",j, sep="")
  }

  dmat = dm.xyz(sim116$trj[i,], grpby=pdb$atom[,"resno"], scut=5)
  contact.apo = dmat.reduce(dmat, 524, 4)
  dmat = dm.xyz(sim129$trj$atp[i,], grpby=pdb$atom[,"resno"], scut=5)
  contact.atp = dmat.reduce(dmat, 524, 4)

  diff=contact.apo-contact.atp
  g3<-which(diff==1, arr.ind=TRUE)
  g4<-which(diff==-1, arr.ind=TRUE)
  
  png(paste("timeDM_", j, ".png", sep=""), width=1000)
  par(mfcol=c(1,2))
  
  g<-which(contact.apo==1, arr.ind=TRUE)
  plot.bio3d(g[,1], g[,2], type='p', cex=cex, pch=pch[1], col=col[1],
             xlim=c(1,524), ylim=c(1,524), ylim2zero=FALSE, 
             sse=sse, sse.border=TRUE, 
             top=top, bot=bot, xlab=xlab, ylab=ylab)

  for (j in 1:nrow(g3))
    draw.circle(g3[j,"row"], g3[j,"col"],  3, border="red", lty="dashed")
  

  g<-which(contact.atp==1, arr.ind=TRUE)
  plot.bio3d(g[,1], g[,2], type='p', cex=cex, pch=pch[2], col=col[2],
             xlim=c(1,524), ylim=c(1,524), ylim2zero=FALSE, 
             sse=sse, sse.border=TRUE, 
             top=top, bot=bot, xlab=xlab, ylab=ylab)

  for (j in 1:nrow(g4))
    draw.circle(g4[j,"row"], g4[j,"col"],  3, border="red", lty="dashed")
  
  dev.off()

}





m1 <- atom.select(pdb, "///33////")
m2 <- atom.select(pdb, "///482////")

m1 <- atom.select(pdb, "///133////")
m2 <- atom.select(pdb, "///474////")

m1 <- atom.select(pdb, "///82////")
m2 <- atom.select(pdb, "///326////")


distance <- function(pdb, trj, sele.1=NULL, sele.2=NULL) {
  
  inds.xyz <- c(sele.1$xyz, sele.2$xyz)

  d <- c()

  for(i in 1:nrow(trj)) {
    dmat <- dist(matrix(trj[i,inds.xyz], ncol = 3, byrow = TRUE))
    dmat = as.matrix(round(dmat,1))
    dmat = dmat[ (length(sele.1$atom)+1) : (length(sele.1$atom)+length(sele.2$atom)),
      1:(length(sele.1$atom))]
    d=c(d, min(dmat))
    
  }
    
  return(d)

}



m1 <- atom.select(pdb, "///82////")
m2 <- atom.select(pdb, "///326////")
m1 <- atom.select(pdb, "///57////")
m2 <- atom.select(pdb, "///208////")


dist.apo=distance(pdb, sim116$trj, m1,m2)
dist.atp=distance(pdb, sim129$trj$atp, m1,m2)


pdf("dist1.pdf", w=12)
par(mfrow=c(3,1), mar=c(2,2,0.5,0.5))
plot(dist.apo, type="o")
lines(dist.atp, col="red")

hist(dist.apo, breaks = 300, main = "", xlim=c(2.5, 8))
hist(dist.atp, breaks = 300, main = "", xlim=c(2.5, 8))
dev.off()



m=matrix(dist.apo, ncol=5)
pdf("mmh.pdf")
par(mfcol=c(2,2))
plot(m[,1], m[,2])##,xlim=c(2,7), ylim=c(2,7))
plot(m[,2], m[,3])##,xlim=c(2,7), ylim=c(2,7))
plot(m[,3], m[,4])##,xlim=c(2,7), ylim=c(2,7))
plot(m[,4], m[,5])##,xlim=c(2,7), ylim=c(2,7))
dev.off()

cor(m[,1], m[,2])







pdf("time_1.pdf")
plot.dcm(dmat, pdb, xlim=c(1,524), ylim=c(1,524),
         sse=sse, sse.grid=TRUE,
         )







md.avg.diff <- sim116$dm$dmat.avg - sim129$dm$atp$dmat.avg


diff <- which((md.avg.diff!=0) & (dmat!=0), arr.ind=T )
#diff <- which((md.avg.diff!=0), arr.ind=T )
diff = cbind(diff, md.avg.diff[diff])
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)),]

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("116_129_dcm.pdf")
plot.dcm(dmat, 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="+")