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

##pdb <- read.pdb("1XCK_reference.pdb")
pdb <- read.pdb("1XCK_chainALM_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 <- "20-50ns_noWAT_300frames_noH_chain_"
trj.inds <- seq(101, 300, by=1)

## inds for statistical tests on distances
##trj.inds <- seq(101, 300, by=10)


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

for ( i in c(1,2,3,4,5,7) ) {
  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_trimer/"
for ( i in c(1,2,3,4,5,7) ) {
  tmptrj <- read.ncdf(paste(prefix129, prefix, i , ".nc", sep=""))
  sim129$trj$atp = rbind(sim129$trj$atp, tmptrj[trj.inds,])
}

sim171b <- NULL
prefix171 <- "/net/lutefisk/slars/groel_md/1XCK_rhodanese/171_1XCK_rhodanese_ATP/results/traj_trimer/"

for ( i in c(1,2,3,4,5,7) ) {
  tmptrj <- read.ncdf(paste(prefix171, prefix, i , ".nc", sep=""))
  sim171b$trj$atp = rbind(sim171b$trj$atp, tmptrj[trj.inds,])
}



sim116$rmsd <- rmsd(pdb$xyz, b=sim116$trj, a.inds=ca.inds$xyz, b.inds=ca.inds$xyz, fit=TRUE)
sim129$rmsd <- rmsd(pdb$xyz, b=sim129$trj$atp, a.inds=ca.inds$xyz, b.inds=ca.inds$xyz, fit=TRUE)
plot(sim116$rmsd)
lines(sim129$rmsd)
dev.off()

sim116$dm <- dm.xyz.trj( sim116$trj, pdb,  threads=6 )
sim129$dm <- dm.xyz.trj( sim129$trj$atp, pdb, threads=6 )


##save(sim116, sim129, file="sim116_129_checkx_no7.RData")
load("sim116_129_checkx_no7.RData")


dmat <- dcm(sim116$dm, sim129$dm, occupancy=0.5, diff.cut=.5)
inds <- which(dmat != 0, arr.ind=T)



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_whole.pdf")
plot.dcm(dmat, pdb, xlim=c(1,1572), ylim=c(1,1572),
         sse=sse, sse.grid=TRUE,
         )
pointLabel(inds[,1], inds[,2], labels, cex=0.3, offset=10)
dev.off()



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

abline(v=524, lty=2)
abline(h=524, lty=2)


dev.off()




md.avg.diff <- sim116$dm$dmat.avg - sim129$dm$dmat.avg
diff <- which((md.avg.diff!=0) & (dmat!=0), arr.ind=T )
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)


t <- which(diff.table1[,"row"]<525 & diff.table1[,"col"]>524)
u <- which(diff.table1[,"col"]<525 & diff.table1[,"row"]>524)
t=c(t,u)



prefix <- "20-50ns_noWAT_300frames_noH_chain_"
trj.inds <- seq(101, 300, by=1)

sim116b <- NULL
prefix116 <- "/net/gulrotkake/slars/groel_md/1XCK/116_1XCK_apo/results/traj_trimer/"

for ( i in c(1,2,3,4,5,7) ) {
  tmptrj <- read.ncdf(paste(prefix116, prefix, i , ".nc", sep=""))
  sim116b$trj = rbind(sim116b$trj, tmptrj[trj.inds,])
}

sim129b <- NULL
prefix129 <- "/net/gulrotkake/slars/groel_md/1XCK/129_1XCK_MGATP/results/traj_trimer/"
for ( i in c(1,2,3,4,5,7) ) {
  tmptrj <- read.ncdf(paste(prefix129, prefix, i , ".nc", sep=""))
  sim129b$trj$atp = rbind(sim129b$trj$atp, tmptrj[trj.inds,])
} 

sim171b <- NULL
prefix171 <- "/net/lutefisk/slars/groel_md/1XCK_rhodanese/171_1XCK_rhodanese_ATP/results/traj_trimer/"

for ( i in c(1,2,3,4,5,7) ) {
  tmptrj <- read.ncdf(paste(prefix171, prefix, i , ".nc", sep=""))
  sim171b$trj$atp = rbind(sim171b$trj$atp, tmptrj[trj.inds,])
}





## check for significance. assume independence
dists.116 <- c()
dists.129 <- c()
dists.171 <- c()

##nrow(diff.table1) ) {
for ( i in 1:nrow(diff.table1) ) {
  mask = paste(diff.table1[i,"row"], diff.table1[i, "col"], sep=",")
  
  inds.tmp <- atom.select(pdb, resno=mask)
  d1 <- dm.inds( sim116b$trj, pdb, inds.tmp)
  d2 <- dm.inds( sim129b$trj$atp, pdb, inds.tmp)
  d3 <- dm.inds( sim171b$trj$atp, pdb, inds.tmp)

  dists.116 <- c(dists.116, d1)
  dists.129 <- c(dists.129, d2)
  dists.171 <- c(dists.171, d3)
}

dists.116.mat <- matrix(dists.116, nrow=nrow(diff.table1), byrow=T)
dists.129.mat <- matrix(dists.129, nrow=nrow(diff.table1), byrow=T)
dists.171.mat <- matrix(dists.171, nrow=nrow(diff.table1), byrow=T)

mean1 <- round(rowMeans(dists.116.mat),1)
mean2 <- round(rowMeans(dists.129.mat),1)
mean3 <- round(rowMeans(dists.171.mat),2)
std1 <- round(apply(dists.116.mat, 1, sd),1)
std2 <- round(apply(dists.129.mat, 1, sd),1)
std3 <- round(apply(dists.171.mat, 1, sd),2)

trj.inds.tmp <- seq(20,1200, by=20) ## every 2nd ns
p.values.t <- apply( cbind(dists.116.mat[,trj.inds.tmp], dists.129.mat[,trj.inds.tmp]), 1,
        function(x) t.test( x[1:(length(x)/2)], x[1+(length(x)/2):length(x)])$p.value )

p.values.w <- apply( cbind(dists.116.mat[,trj.inds.tmp], dists.129.mat[,trj.inds.tmp]), 1,
        function(x) wilcox.test( x[1:(length(x)/2)], x[1+(length(x)/2):length(x)])$p.value )
p.values.w.3 <- apply( cbind(dists.116.mat[,trj.inds.tmp], dists.171.mat[,trj.inds.tmp]), 1,
        function(x) wilcox.test( x[1:(length(x)/2)], x[1+(length(x)/2):length(x)])$p.value )


p.values.t <- round(p.values.t, 5)
p.values.w <- round(p.values.w, 5)
p.values.w.3 <- round(p.values.w.3, 5)

results <- cbind(diff.table1, mean1, mean2, mean3, std1, std2, std3, p.values.w, p.values.w.3)


x <- which(results[,"row"]<525)
y <- which(results[,"col"]>524)
contact.inds <- intersect(x,y)

sign.inds <- which(p.values.w<0.05)
sign.inds <- intersect(sign.inds, which(p.values.w.3<0.05))
results[sign.inds,]

write.table(results, file="statiscs_distances.dat")
save(results, sign.inds, file="final_results_statistics.RData")



## HOLO contacts
holo.inds <- which(results[,"avg.diff"]>0)
holo.inds <- intersect(holo.inds, sign.inds)
tmp1=results[holo.inds,]
c1 <- paste(tmp1[,"row"], collapse="+")
c2 <- paste(tmp1[,"col"], collapse="+") 
paste(c1,c2,sep="+")
write.table(tmp1[,c("avg.diff", "res1", "res2")], file="holo_contacts.dat")
  
## APO contacts
apo.inds <- which(results[,"avg.diff"]<0)
apo.inds <- intersect(apo.inds, sign.inds)
tmp2=results[apo.inds,]
c1 <- paste(tmp2[,"row"], collapse="+")
c2 <- paste(tmp2[,"col"], collapse="+")
paste(c1,c2,sep="+")
write.table(tmp2[,c("avg.diff", "res1", "res2")], file="apo_contacts.dat")





hmm <- c('from pymol.cgo import *')
for ( i in 1:nrow(tmp1) ) {
  name = paste('contact', i, sep="")
  h = write.interaction(tmp1[i,"row"], tmp1[i,"col"], pdb, name, color=c(1,0,0))
  hmm=c(hmm,h)
}

for ( j in 1:nrow(tmp2) ) {
  name = paste('contact', i+j, sep="")
  h = write.interaction(tmp2[j,"row"], tmp2[j,"col"], pdb, name, color=c(0,0,1))
  hmm=c(hmm,h)
}
write(hmm, file="test.py")


write.interaction <- function(a,b,pdb, name, color) {

  resa <- atom.select(pdb, resno=a, elety="CA")
  resb <- atom.select(pdb, resno=b, elety="CA")
  
  x1 = pdb$atom[resa$atom,"x"]
  y1 = pdb$atom[resa$atom,"y"]
  z1 = pdb$atom[resa$atom,"z"]

  x2 = pdb$atom[resb$atom,"x"]
  y2 = pdb$atom[resb$atom,"y"]
  z2 = pdb$atom[resb$atom,"z"]

  r1 = color[1]
  g1 = color[2]
  b1 = color[3]
  
  radius = 0.4

  name = paste("'", name, "'", sep="")
  nums =  paste( '9.0', x1, y1, z1, x2, y2, z2, radius, r1, g1, b1, r1, g1, b1, sep=", ")
  str = paste("cmd.load_cgo( [ ", nums, "], ", name, " )", sep=" " )

  return(str)
  
}