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

##pdb <- read.pdb("1XCK_reference.pdb")
pdb <- read.pdb("1XCK_apo_dimer_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 <- "30-50ns_noWAT_noH_1000frames_chain_"
trj.inds <- seq(1, 1000, by=100)


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

for ( i in 1:6 ) {
  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_dimers/"


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


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


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


dmat <- dcm(sim116$dm, sim129$dm$atp, occupancy=0.5, diff.cut=0.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.pdf")
plot.dcm(dmat, pdb, xlim=c(1,1048), ylim=c(1,1048),
         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=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$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)



prefix <- "30-50ns_noWAT_noH_1000frames_chain_"
trj.inds <- seq(5, 1000, by=5)


sim116a <- NULL
prefix116 <- "/net/gulrotkake/slars/groel_md/1XCK/116_1XCK_apo/results/traj_dimers/"

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

sim116b <- NULL
prefix116 <- "/net/gulrotkake/slars/groel_md/1XCK/116_1XCK_apo/results/traj_dimers/"
prefix <- "20-50ns_noWAT_1500frames_noH_chain_"
trj.inds <- seq(505, 1500, by=5)

for ( i in 8:13 ) {
  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_dimers/"
prefix <- "30-50ns_noWAT_noH_1000frames_chain_"
trj.inds <- seq(5, 1000, by=5)

for ( i in 1:6 ) {
  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_dimers/"
prefix <- "20-50ns_noWAT_1500frames_noATP_noH_chain_"
trj.inds <- seq(505, 1500, by=5)

for ( i in 1:6 ) {
  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.116b <- c()
dists.129 <- c()
dists.171 <- c()
dists.apo <- c()
dists.atp <- c()
dists.adp <- c()

for ( i in 1:nrow(diff.table1) ) {
  mask = paste(diff.table1[i,"row"], diff.table1[i, "col"], sep=",")

  if ( diff.table1[i, "avg.diff"] > 0 )
    alternative = "greater"
  else
    alternative = "less"
  
  inds.tmp <- atom.select(pdb, resno=mask)
  d1 <- dm.inds( sim116a$trj, pdb, inds.tmp)
  d1b <- 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.116b <- c(dists.116b, d1b)
  dists.129 <- c(dists.129, d2)
  dists.171 <- c(dists.171, d3)

  #w <- wilcox.test(d1,d2, paired = FALSE, alternative=alternative)
  #sign = c(sign, w$p.value)
}

dists.116.mat <- matrix(dists.116, nrow=nrow(diff.table1), byrow=T)
dists.116b.mat <- matrix(dists.116b, 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),2)
mean1b <- round(rowMeans(dists.116b.mat),2)
mean2 <- round(rowMeans(dists.129.mat),2)
mean3 <- round(rowMeans(dists.171.mat),2)

std1 <- round(apply(dists.116.mat, 1, sd),2)
std1b <- round(apply(dists.116b.mat, 1, sd),2)
std2 <- round(apply(dists.129.mat, 1, sd),2)
std3 <- round(apply(dists.171.mat, 1, sd),2)


# dont use all frames for statistical testing
#trj.inds.tmp <- seq(10,1400, by=10) ## every 1nd ns
trj.inds.tmp <- seq(20,1200, by=20) ## every 2nd ns
#trj.inds.tmp <- seq(50,1400, by=50) ## every 5th ns

# indices for every second subunit, using every 2nd ns
x <- seq(20,200, by=20)
trj.inds.tmp.a <- c(x, x+400, x+800)
trj.inds.tmp.b <- c(x+200, x+600, x+1000)
trj.inds.tmp.c <- c(x+400, x+800, x+1200)


# right? 
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.116b.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, mean1b, mean2, mean3, std1, std1b, std2, std3, p.values.w, p.values.w.3)
results[,6:14]

sign.inds <- which(p.values.w<0.05)
sign.inds <- intersect(sign.inds, which(p.values.w.3<0.05))
results[sign.inds,c("avg.diff", "res1", "res2", "mean1", "std1",  "mean2", "std2", "mean3", "p.values.w.3")]

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


results[sign.inds,c("avg.diff", "res1", "res2", "mean1", "mean1b", "mean2", "mean3", "dists.apo", "dists.atp", "dists.adp", "p.values.w", "p.values.w.3")]
results[sign.inds,c("avg.diff", "res1", "res2", "mean1", "std1", "mean1b", "std1b", "mean2", "std2", "mean3", "std3")]

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





## HOLO contacts
#tmp1=tail(diff.table1, n=50)
holo.inds <- which(diff.table1[,"avg.diff"]>0)
holo.inds = intersect(holo.inds, sign.inds)
tmp1=diff.table1[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
#tmp2=head(diff.table1, n=50)
apo.inds <- which(diff.table1[,"avg.diff"]<0)
apo.inds = intersect(apo.inds, sign.inds)
tmp2=diff.table1[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)
  
}