source("mmpbsa_funs.R")
library(bio3d)
library(Hmisc)
pdb <- read.pdb("1XCK_chainA_noWAT_noH.pdb")
seq <- seq.pdb(pdb)
s <- array(seq)



a1 <- read.mmpbsa.decomp("116/concat_FINAL_DECOMP_MMPBSA.dat")
b1 <- read.mmpbsa.decomp("129/concat_FINAL_DECOMP_MMPBSA.dat")
c1 <- read.mmpbsa.decomp("171/FINAL_DECOMP_MMPBSA_1.dat")
#d <- read.mmpbsa.decomp("191/FINAL_DECOMP_MMPBSA_1.dat")

apo <- calc.means(list(a1), n=1)
holo <- calc.means(list(b1), n=1)
rhod <- calc.means(list(c1), n=1)
#e434k <- calc.means(list(d), n=1)


# this is for importing to PYMOL
inds <- 1:1572
s.inds <- ((c(1:1572)-1)%%524)+1
hmm <- cbind(paste(s[s.inds], inds, sep=""), apo$total)
write.table(hmm, file="apo.total.dat", row.names=F, col.names=F)



#here goes the plotting of residue contributions
inds <- which(abs(apo$total)>0.5)
#inds=inds[which(inds<525)]
cex=1.4
cex.legend=1.0
cex.axis=1.4
cex.label=1.4

n <- 3600
ylim=c(-6.5,3)
#diff.inds <- which(abs(apo$total-holo$total)>0.5)
#diff.inds=diff.inds[1:6]

# use these for WT (1XCK)
diff.inds=c(inds[1:4],433,inds[5:length(inds)])
#diff.inds=inds  ## unchanged
seq <- seq.pdb(pdb)
grps=c(10,5,5)
helplines=seq(-5,1)

if ( FALSE ) {
# use these for E434K
diff.inds=inds
diff.inds=c(inds[1:2], 429, inds[3:17],1477,inds[18:length(inds)])
seq <- seq.pdb(pdb)
seq[433]="K"
grps=c(10,5,6)
}

filename="residues_with_rhodCheck.pdf"

grps=cumsum(grps)

pdf(filename, w=8, h=10)
par(mfcol=c(2,1), mar=c(3,3,2,1), mgp=c(1.75,.5,0))
mp <- plot.zoom.three(apo, holo, rhod, diff.inds, n, seq, ylim, 
                mtext.line=-0.25, cex=cex, cex.axis=cex.axis, cex.lab=cex.label)
mtext(3, at=0, text="A", cex=2.5)

##abline(h=helplines, lty=3, col="grey90")

a=1;b=grps[1];
y=c(ylim[1], ylim[1])+0.5
mtext(1, at=mean(c(mp[1,a], mp[1,b])), adj=0.5, line=-2.2, cex=cex.axis, text="cis subunit")
lines(x=c(mp[1,a], mp[2,b]), y=y)
a=grps[1]+1;b=grps[2];
mtext(1, at=mean(c(mp[1,a], mp[1,b])), adj=0.5, line=-2.2, cex=cex.axis, text="trans (right)")
lines(x=c(mp[1,a], mp[2,b]), y=y)
a=grps[2]+1;b=grps[3];
mtext(1, at=mean(c(mp[1,a], mp[1,b])), adj=0.5, line=-2.2, cex=cex.axis, text="trans (left)")
lines(x=c(mp[1,a], mp[2,b]), y=y)

legend("topright", legend=c("GroEL",
                  expression(paste("GroEL", " ATP"[7])),
                  expression(paste("GroEL"[Rhod], "ATP"[7]))),
       fill=c("gray20", "gray50", "gray90"), cex=cex.legend, ncol=3)



#cex=.8
ylim=ylim
s.inds <- ((diff.inds-1)%%524)+1
hmm <- NULL
hmm = t(cbind(apo$vdw,  holo$vdw, rhod$vdw, apo$ele, holo$ele, rhod$ele, apo$sas, holo$sas, rhod$sas))

mp <- barplot(hmm[,diff.inds], beside=T, ylim=ylim,
              ylab="Energy (kcal/mol)", col=c(rep("gray20",3),rep("gray50",3), rep("gray90",3) ),
              cex=cex, cex.axis=cex.axis, cex.lab=cex.label)
mtext(1, line=-0.25, at=colMeans(mp), las=2, cex=cex,
      text=paste(seq[s.inds], s.inds+1))
mtext(3, at=0, text="B", cex=2.5)

a=1;b=grps[1];
y=c(ylim[1],ylim[1])+0.5
mtext(1, at=mean(c(mp[1,a], mp[nrow(mp),b])), adj=0.5, line=-2.2, cex=cex.axis, text="cis subunit")
lines(x=c(mp[1,a], mp[nrow(mp),b]), y=y)
a=grps[1]+1;b=grps[2];
mtext(1, at=mean(c(mp[1,a], mp[nrow(mp),b])), adj=0.5, line=-2.2, cex=cex.axis, text="trans (right)")
lines(x=c(mp[1,a], mp[nrow(mp),b]), y=y)
a=grps[2]+1;b=grps[3];
mtext(1, at=mean(c(mp[1,a], mp[nrow(mp),b])), adj=0.5, line=-2.2, cex=cex.axis, text="trans (left)")
lines(x=c(mp[1,a], mp[nrow(mp),b]), y=y)

legend("topright", legend=c("VDW", "ELE", "SAS"), cex=cex.legend,
       fill=c("gray20", "gray50", "gray90"), ncol=3)
dev.off()
















## the rest here is just testing
## and having fun

if(FALSE){
l.inds <- which(abs(a1$deltas$total)>1.5)
l.inds2 <- ((l.inds-1)%%524)+1

par(mfcol=c(2,1))
mp <- barplot(a1$deltas$total)
text(mp[l.inds], a1$deltas$total[l.inds], labels=paste(s[l.inds2], l.inds, sep=""), las=2)

l.inds <- which(abs(b1$deltas$total)>1.5)
l.inds2 <- ((l.inds-1)%%524)+1


mp <- barplot(b1$deltas$total)
text(mp[l.inds], b1$deltas$total[l.inds], labels=paste(s[l.inds2], l.inds, sep=""), las=2)
}



ifr <- c(21,24:25,28,33:50,55:63,456:460,480:482, 382:388, 177:182, 215:216, 267:271, 240:245, 320, 641, 635:637, 1036:1048, 525:539, 588:603, 751:752, 754:755, 778:780, 748:749, 800:804, 719:720, 807, 881, 883, 887)
ifr.rec <- c(21,24:25,28,33:50,55:63,456:460,480:482, 382:388, 177:182, 215:216, 267:271, 240:245, 320)
ifr.lig <- c(641, 635:637, 1036:1048, 525:539, 588:603, 751:752, 754:755, 778:780, 748:749, 800:804, 719:720, 807, 881, 883, 887)
#ifr.lig=((ifr.lig-1)%%524)+1


cutoff <- 1
apo.means <- colMeans(rbind(a1$deltas$total, a2$deltas$total, a3$deltas$total))
holo.means <- colMeans(rbind(b1$deltas$total, b2$deltas$total, b3$deltas$total))
diff <- apo.means-holo.means



ylim=c(-8,3.4)
offset=.5
label.cex=.8
axis.cex=.8
ylab="Energy (kcal/mol)"

pdf("means_mmpbsa_decom.pdf", w=12, h=6)
par(mfrow=c(4,1), mar=c(1,3,1,1), mgp=c(1.75,0.75,0))

l.inds <- which(apo.means[1:524]>1.5)
s.inds <- which(holo.means[1:524]<(-1.5))
l.inds=c(33,50,74,215,389)
s.inds=c(38,46,182,269,385)

mp <- barplot(apo.means[1:524], ylim=ylim, ylab=ylab)
mtext(3, at=-20, adj=0, text="A", cex=1.75, line=-1)
text(mp[s.inds], apo.means[s.inds]-offset, labels=paste(s[s.inds], s.inds+1, sep=""), las=2, cex=label.cex)
text(mp[l.inds], apo.means[l.inds]+offset, labels=paste(s[l.inds], l.inds+1, sep=""), las=2, cex=label.cex)


l.inds=c(l.inds, 41,244)
s.inds=c(38,46,182,241,267,385)


mp <- barplot(holo.means[1:524], ylim=ylim, ylab=ylab)
mtext(3, at=-20, adj=0, text="B", cex=1.75, line=-1)
text(mp[s.inds], holo.means[s.inds]-offset, labels=paste(s[s.inds], s.inds+1, sep=""), las=2, cex=label.cex)
text(mp[l.inds], holo.means[l.inds]+offset, labels=paste(s[l.inds], l.inds+1, sep=""), las=2, cex=label.cex)

dev.off()



#pdf("lig_means_mmpbsa_decom.pdf", w=12, h=6)
#par(mfrow=c(2,1), mar=c(1,3,1,1), mgp=c(1.75,0.75,0))
inds <- 525:1048

l.inds <- which(apo.means[inds]>1)
s.inds <- which(holo.means[inds]<(-1))
l.inds=c(12,116,224,251,302)
s.inds=c(7,68,113,196,230,280,284,359,518)

mp <- barplot(apo.means[inds], ylim=ylim, ylab=ylab)
mtext(3, at=-20, adj=0, text="A", cex=1.75, line=-1)
text(mp[s.inds], apo.means[s.inds+524]-offset, labels=paste(s[s.inds], s.inds+1, sep=""), las=2, cex=label.cex)
text(mp[l.inds], apo.means[l.inds+524]+offset, labels=paste(s[l.inds], l.inds+1, sep=""), las=2, cex=label.cex)


#l.inds=c(33,50,74,215,389)
#s.inds=c(38,182,241, 267,385)

mp <- barplot(holo.means[inds], ylim=ylim, ylab=ylab)
mtext(3, at=-20, adj=0, text="B", cex=1.75, line=-1)
text(mp[s.inds], holo.means[s.inds+524]-offset, labels=paste(s[s.inds], s.inds+1, sep=""), las=2, cex=label.cex)
text(mp[l.inds], holo.means[l.inds+524]+offset, labels=paste(s[l.inds], l.inds+1, sep=""), las=2, cex=label.cex)

dev.off()


























pdf("diff_mmpbsa_decomo.pdf", w=12, h=6)
par(mfcol=c(3,1), mar=c(3,2,1,1))
plot(apo.means, type="o", cex=0.5)
points(ifr, apo.means[ifr], col="green", cex=0.5)

plot(holo.means, type="o", cex=0.5)
points(ifr, holo.means[ifr], col="green", cex=0.5)

plot(diff, type="o", cex=0.5)
points(ifr, diff[ifr], col="green", cex=0.5)
peaks=which(abs(diff)>cutoff)
p.inds=((peaks-1)%%524)+1
text(peaks, diff[peaks]*1.05, labels=paste(s[p.inds], peaks), col="red")
dev.off()








par(mfcol=c(3,1))
plot(a1$deltas$total, type="o", cex=0.5)
peaks=which(abs(a1$deltas$total)>2)
p.inds=((peaks-1)%%524)+1
text(peaks, a1$deltas$total[peaks]*1.05, labels=paste(s[p.inds], peaks), col="red")
points(ifr, a1$deltas$total[ifr], col="green", cex=0.5)


plot(b1$deltas$total, type="o", cex=0.5)
peaks=which(abs(b1$deltas$total)>2)
p.inds=((peaks-1)%%524)+1
text(peaks, b1$deltas$total[peaks]*1.05, labels=paste(s[p.inds], peaks), col="red")
points(ifr, b1$deltas$total[ifr], col="green", cex=0.5)

plot(c1$deltas$total, type="o", cex=0.5)
peaks=which(abs(c1$deltas$total)>2)
p.inds=((peaks-1)%%524)+1
text(peaks, c1$deltas$total[peaks]*1.05, labels=paste(s[p.inds], peaks), col="red")
points(ifr, c1$deltas$total[ifr], col="green", cex=0.5)