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

##pdb <- read.pdb("1XCK_reference.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_"

# for contact map calculation, use frames every 100 ps
trj.inds <- seq(3010, 5000, by=10)

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

for ( i in 1: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_monomer/"


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


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

sim116$trj <- NULL
sim129$trj$atp <- NULL
#save(sim116, sim129, file="sim116_129_new.RData")
load("sim116_129_new.RData")

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

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


diff.inds <- which((md.avg.diff!=0) & (dmat!=0), arr.ind=T )
#diff <- which((md.avg.diff!=0), arr.ind=T )
diff = cbind(diff.inds, md.avg.diff[diff.inds])
diff = cbind(diff, sim116$dm$dmat.avg[diff.inds])
diff = cbind(diff, sim129$dm$atp$dmat.avg[diff.inds])

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", "d1", "d2", "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="+")




## using all 5000 frames - noH
prefix <- "50ns_noWAT_5000frames_noH_chain_"

## for test of significane use only last 20 ns, 2 ns interval
#trj.inds <- seq(3200, 5000, by=200)

## to report distances, use 0.2 ns interval
trj.inds <- seq(3010, 5000, by=10)


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

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

for ( i in 8:14 ) {
  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_monomer/"


for ( i in 1: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_monomer/"

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

pdb.apo <- read.pdb("1XCK_A.pdb")
pdb.atp <- read.pdb("1SX3_A.pdb")
pdb.adp <- read.pdb("1SX4_A.pdb")
dm.apo <- dm.xyz(pdb.apo$xyz, grpby=pdb.apo$atom[,"resno"])
dm.atp <- dm.xyz(pdb.atp$xyz, grpby=pdb.atp$atom[,"resno"])
dm.adp <- dm.xyz(pdb.adp$xyz, grpby=pdb.adp$atom[,"resno"])

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)
  #trj.inds.tmp <- seq(20,1400, by=20)
  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)

  dists.apo <- c(dists.apo, dm.apo[diff.table1[i,"row"],diff.table1[i, "col"]])
  dists.atp <- c(dists.atp, dm.atp[diff.table1[i,"row"],diff.table1[i, "col"]])
  dists.adp <- c(dists.adp, dm.adp[diff.table1[i,"row"],diff.table1[i, "col"]])

  #t <- t.test(d1[trj.inds.tmp],d2[trj.inds.tmp])
  #sign.t <- c(sign.t, t$p.value)

  #w <- wilcox.test(d1[trj.inds.tmp],d2[trj.inds.tmp], 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)

dists.apo <- round(dists.apo, 2)
dists.atp <- round(dists.atp, 2)
dists.adp <- round(dists.adp, 2)


# dont use all frames for statistical testing
#trj.inds.tmp <- seq(10,1400, by=10) ## every 1nd ns
trj.inds.tmp <- seq(20,1400, 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.w.a <- apply( cbind(dists.116.mat[,trj.inds.tmp.a], dists.129.mat[,trj.inds.tmp.a]), 1,
        function(x) wilcox.test( x[1:(length(x)/2)], x[1+(length(x)/2):length(x)])$p.value )
p.values.w.b <- apply( cbind(dists.116.mat[,trj.inds.tmp.b], dists.129.mat[,trj.inds.tmp.b]), 1,
        function(x) wilcox.test( x[1:(length(x)/2)], x[1+(length(x)/2):length(x)])$p.value )
p.values.w.c <- apply( cbind(dists.116.mat[,trj.inds.tmp.c], dists.129.mat[,trj.inds.tmp.c]), 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)

p.values.w.a <- round(p.values.w.a, 5)
p.values.w.b <- round(p.values.w.b, 5)
p.values.w.c <- round(p.values.w.c, 5)

results <- cbind(diff.table1, mean1, mean1b, mean2, mean3, std1, std1b, std2, std3, dists.apo, dists.atp, dists.adp, p.values.w, p.values.w.3, p.values.w.a, p.values.w.b, p.values.w.c) 
results[,6:14]

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

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
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)
  
}