library(bio3d)
library(Hmisc)


##pdb <- read.pdb("/net/lutefisk/slars/groel_md/bio3d_groel/raw_pdbs/split_chain/1XCK_A.pdb")
pdb <- read.pdb("~/lutefisk/groel_md/bio3d_groel/raw_pdbs/split_chain/1XCK_A.pdb")
sse <- dssp(pdb)


vlines=c(132,408,190,376)
hlines=c(132,408,190,376)

corr.MD.holo <- as.matrix(read.table("/home/slars/gulrotkake/groel_md/1XCK/129_1XCK_MGATP/results/correlation_dimers/dimer_corr.dat"))
corr.MD.apo<- as.matrix(read.table("/home/slars/gulrotkake/groel_md/1XCK/116_1XCK_apo/results/correlation_dimers/dimer_corr.dat"))











pdf("corrMat.holo.pdf", w=7.6, h=7)
plot.corrmat(corr.MD.holo[525:1048,1:524], sse, vlines, hlines, main="Cross correlation")
dev.off()

pdf("corrMat.apo.pdf", w=7.6, h=7)
plot.corrmat(corr.MD.holo[525:1048,1:524], sse, vlines, hlines, main="Cross correlation")
dev.off()









plot.corrmat <- function(mat, sse=NULL, vlines=NULL, hlines=NULL,
                         xlab = "Residue index",
                         ylab = "Residue index",
                         main = "Correlation plot",
                         levels = pretty( c(-1,1), 50),
                         Labpalette = colorRampPalette(c("blue", "white", "red"),
                           space = "Lab"),
                         cex.axis=.7, cex.sse=.6, ...) {


  nr <- dim(mat)[1]
  nc <- dim(mat)[2]



  filled.contour( x=1:nr, y=1:nc, z=mat, color = Labpalette,
                 levels = levels,
                 nlevels = 10,
                 xlab = xlab, ylab=ylab, main = main,

                 plot.axes={
                   ##axis(1); axis(2);
                   axis(2, at=sse$helix$start, labels=sse$helix$start, cex.axis=cex.axis);
                   axis(1, at=sse$helix$start, labels=sse$helix$start, cex.axis=cex.axis);
                   abline(v=sse$helix$start, lty=2, col="grey70");
                   abline(h=sse$helix$start, lty=2, col="grey70");
                   mtext(capitalize(letters[1:length(sse$helix$start)]), side=4, line=0.2, las=2, at=(sse$helix$start+sse$helix$end)/2, cex=cex.sse);
                   mtext(capitalize(letters[1:length(sse$helix$start)]), side=3, line=0.2, at=(sse$helix$start+sse$helix$end)/2, cex=cex.sse);
                   abline(v=vlines); abline(h=hlines);
                   mtext(c("E1", "I1", "A", "I2", "E2"), at=c(132/2, (132+190)/2, (190+376)/2, (376+408)/2, (408+524)/2), side=3, line=0.7, cex=1.1);

                   abline(0,1);

                   #rect(25,376,55,408, border="green", lwd=3);
                   #rect(5,325,90,345, border="green", lwd=3);
                   #rect(136,340,190,376, border="green", lwd=3);
                   #rect(340,376,376,408, border="green", lwd=3);
                   #rect(25,135,55,190, border="green", lwd=3);

                   #rect(136,25,190,60, border="yellow", lwd=3);
                   #rect(376,25,408,60, border="yellow", lwd=3);
                   #rect(462,340,500,376, border="yellow", lwd=3);

                 },
                 ...
                 )
}