library(bio3d)
library(Hmisc)
source("corr_mat_fun.R")

pdb.A <- read.pdb("1XCK_A.pdb")
sse.A <- dssp(pdb.A)

pdb <- read.pdb("trimer_ADP_2xATP.pdb")
sse <- dssp(pdb)

#vlines=c(525,1048,1572,1584,1597)
#hlines=c(525,1048,1572,1584,1597)


domains=c(137,410,375,190)
domains2=c(domains+524)
vlines=c(525,1048,1060,1073, domains, domains2)
hlines=c(domains,525,536,549)


helix.labels.y = capitalize(letters[1:length(sse.A$helix$start)])
helix.labels.x = rep(helix.labels.y, 2)
#helix.labels.x = c(helix.labels.x, "cis ATP", "trans ADP", "trans ADP")

helix.labels.at.y=(sse.A$helix$start+sse.A$helix$end)/2
helix.labels.at.x=(sse$helix$start[1:36]+sse$helix$end[1:36])/2
#helix.labels.at.x=c(helix.labels.at.x, (1048+6), (1048+19), (1048+28))

axis.at.y=sse.A$helix$start
axis.labels.y=sse.A$helix$start

axis.at.x=sse$helix$start[1:36]
axis.labels.x=rep(sse.A$helix$start,2)


#x=1:1601
#y=1:1610
#x=525:1572
x=525:1610
#x=525:1585
#y=1:524
y=c(1:524,1573:1610)

dim=c(length(x),length(y))
dat <- array(0, dim=c(dim[1], dim[2], 12))


apo <- matrix(0,dim[1],dim[2])
j=1
for ( i in c(2,3,4,5,6,7) ) {
  A <- scan(paste("/net/gulrotkake/slars/groel_md/1SX4/180_1SX4/results/corr_trimer/corr_50ns_CA_chain_",i,sep=""))
  n <- sqrt(length(A))
  A <- matrix( A, n, n, byrow = TRUE)
  dat[, ,j]=A[x,y]
  apo=apo+A[x,y]
  j=j+1
}
apo=round(apo/6,2)
##write.table(apo, file="corrmat_apo_avg.dat", row.names=F, col.names=F, quote=F)

##x=1:1610
##y=1:1610
#x=525:1610
#y=1:524

## atp:1049-1060, adp: 1061-1073, 1074-1086
## atp:525-536, adp: 537-549,550-562
## correlation between ATP and ADP:
## holo[1061:1073, 525:536]

c <- holo[1061:1073, 525:536]
pdf("nucleotide_correlation.pdf", w=10, h=4)
par(mfcol=c(2,4), mar=c(1.,1,0.5,0.5))
image(1:13, 1:12, c, col=gray(10:0/10), zlim=c(0,.5),
      xlab="MD-PCs", ylab="X-ray PCs")
#zlim=c(0,1), 
dev.off()





holo <- matrix(0,dim[1],dim[2])
j=1
for ( i in c(2,3,4,5,6,7) ) {
  A <- scan(paste("/net/gulrotkake/slars/groel_md/1SX4/181_1SX4_atp/results/corr_trimer/corr_50ns_CA_chain_",i,sep=""))
  n <- sqrt(length(A))
  print(n)
  A <- matrix( A, n, n, byrow = TRUE)
  dat[, ,j]=A[x,y]
  holo=holo+A[x,y]
  j=j+1
}
holo=round(holo/6,2)

##write.table(holo, file="corrmat_holo_avg.dat", row.names=F, col.names=F, quote=F)

diff=apo-holo
#diff=filter.corrmat(diff, threshold=-1)

v <- t.list
##v[upper.tri(v, diag=FALSE)] <- NA4
t <- which(v<0.005, arr.ind=T)
sig.mat <- matrix(0,dim[1],dim[2])
sig.mat[t]=diff[t]
g <- which(sig.mat!=0, arr.ind=T)



png(paste("diff.png",sep=""), width=2000, height=1000,
    pointsize=24)
plot.corrmat(sig.mat, sse,
##plot.corrmat(diff, sse,             
             helix.labels.x=helix.labels.x, helix.labels.at.x=helix.labels.at.x,
             helix.labels.y=helix.labels.y, helix.labels.at.y=helix.labels.at.y,
             axis.at.x=axis.at.x, axis.labels.x=axis.labels.x,
             axis.at.y=axis.at.y, axis.labels.y=axis.labels.y,
             vlines=vlines, hlines=hlines, main="")
dev.off()

source("corr_mat_fun.R")
A=filter.corrmat(apo, threshold=.0)
png(paste("apo_avg.png",sep=""), width=2000, height=1000,
    pointsize=24)
plot.corrmat(A, sse,
             helix.labels.x=helix.labels.x, helix.labels.at.x=helix.labels.at.x,
             helix.labels.y=helix.labels.y, helix.labels.at.y=helix.labels.at.y,
             axis.at.x=axis.at.x, axis.labels.x=axis.labels.x,
             axis.at.y=axis.at.y, axis.labels.y=axis.labels.y,
             vlines=vlines, hlines=hlines, main="")
dev.off()

source("corr_mat_fun.R")
B=filter.corrmat(holo, threshold=.0)
png(paste("holo_avg.png",sep=""), width=2000, height=1000,
    pointsize=12, res=200)
#pdf("holo_avg.pdf", w=16, h=8)
plot.corrmat(B, sse,
             helix.labels.x=helix.labels.x, helix.labels.at.x=helix.labels.at.x,
             helix.labels.y=helix.labels.y, helix.labels.at.y=helix.labels.at.y,
             axis.at.x=axis.at.x, axis.labels.x=axis.labels.x,
             axis.at.y=axis.at.y, axis.labels.y=axis.labels.y,
             vlines=vlines, hlines=hlines, main="")
dev.off()

png(paste("corrMat.apo_",i,".png",sep=""), width=1500, height=1500,
      pointsize=24)
  plot.corrmat(A, sse,
               helix.labels=helix.labels, helix.labels.at=helix.labels.at,
               axis.at=axis.at, axis.labels=axis.labels,
               vlines=vlines, hlines=hlines, main="Cross correlation")
dev.off()

png(paste("corrMat.holo_",i,".png",sep=""), width=1500, height=1500,
    pointsize=24)
plot.corrmat(A, sse,
             helix.labels=helix.labels, helix.labels.at=helix.labels.at,
             axis.at=axis.at, axis.labels=axis.labels,
             vlines=vlines, hlines=hlines, main="Cross correlation")
  dev.off()




library(multicore)
threads=4
ptm <- proc.time()
rnames=rep(1:threads, each=(nrow(dat)/threads))
print(rnames)

jobs <- list()
for ( k in 1:threads ) {
  r.inds=which(rnames==k)
  mat=dat[r.inds,,]
  q=parallel(t.test.mat(mat))
  jobs[[k]]=q
}

res=collect(jobs, wait=TRUE)
print(proc.time() - ptm)

t.list <- c()
for ( job in res ) {
  t.list=rbind(t.list, job)
}


t.test.mat <- function(mat) {
  n <- dim(mat)[1]
  m <- dim(mat)[2]
  t.list <- matrix(0,n,m)
  
  for ( i in 1:n ) {
    for ( j in 1:m ) {
      a <- mat[i,j,1:6]
      b <- mat[i,j,7:12]

      if (!identical(a,b)) {
        t <- t.test(a,b)
        t.list[i,j]=t$p.value
      }
      else
        t.list[i,j]=1
    }
  }
  return(t.list)
}


filter.corrmat <- function(A, threshold=0.2) {

  h <- which(A>threshold)
  l <- which(A<(threshold*(-1)))
  hmm=A*0
  hmm[h]=A[h]
  hmm[l]=A[l]
  return(hmm)

}