library(bio3d)
library(Hmisc)
source("corr_mat_fun.R")
##source("~/Documents/bio3d_new_funs/plot.cij.R")

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

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


#A <- scan( "/net/gulrotkake/slars/groel_md/1XCK/116_1XCK_apo/results/correlation_trimer/trimer_corr.dat" )
A <- scan( "~/gulrotkake/groel_md/1XCK/116_1XCK_apo/results/correlation_trimer/trimer_corr_all.dat" )
n <- sqrt(length(A))
corr.MD.apo <- matrix( A, n, n, byrow = TRUE)

##A <- scan( "/net/gulrotkake/slars/groel_md/1XCK/129_1XCK_MGATP/results/correlation_trimer/trimer_corr_wATP.dat")
A <- scan( "~/gulrotkake/groel_md/1XCK/129_1XCK_MGATP/results/correlation_trimer/trimer_corr_all.dat")
n <- sqrt(length(A))
corr.MD.holo <- matrix( A, n, n, byrow = TRUE)

h <- which(corr.MD.holo>0.2)
l <- which(corr.MD.holo<(-0.2))
hmm=corr.MD.holo*0
hmm[h]=corr.MD.holo[h]
hmm[l]=corr.MD.holo[l]
corr.MD.holo=hmm


vlines=c(525,1048)
hlines=c(525,1048)


helix.labels = capitalize(letters[1:length(sse.A$helix$start)])
helix.labels=rep(helix.labels, 3)

helix.labels.at=(sse$helix$start+sse$helix$end)/2
axis.at=sse$helix$start[1:54]
axis.labels=rep(sse.A$helix$start,3)
  
##pdf("correlation_plotATP.pdf", w=15.2, h=7)


png("corrMat.holo.png", width=2000, height=1000,
    pointsize=96)
plot.corrmat(corr.MD.holo[525:1585,1:524], 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("corrMat.holo.png", width=2000, height=2000,
    pointsize=24)
plot.corrmat(corr.MD.holo, 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("corrMat.apo.png", width=2000, height=2000,
    pointsize=24)
plot.corrmat(corr.MD.apo, 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()