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

#dm <- dm.xyz(pdb.dimer$xyz, grpby=pdb.dimer$atom[,"resno"], scut=NULL)
#dm=dm[1:524,525:1048]
#ifr <- which(dm < 5, arr.ind=TRUE)


a1 <- read.mmpbsa.decomp("116/FINAL_DECOMP_MMPBSA_1.dat")
b1 <- read.mmpbsa.decomp("129/FINAL_DECOMP_MMPBSA_1.dat")


apo <- calc.means(list(a1), n=1)
holo <- calc.means(list(b1), n=1)





cex=.7
cex.axis=.8
cex.label=.8
n <- 200
ylim=c(-8,3)
diff.inds <- which(abs(apo$total-holo$total)>1.5)


pdf("decomp_desidues.pdf", w=12, h=3)
par(mfcol=c(1,2), mar=c(3,3,2,1), mgp=c(1.75,.6,0))
mp <- plot.zoom(apo, holo, diff.inds, n, seq, ylim,
                cex=cex, cex.axis=cex.axis, cex.lab=cex.label)
mtext(3, at=0, text="A", cex=1.1)
#mtext(1, line=3, at=c(mean(mp[,1:8]), mean(mp[,9:15])),
#      text=c("Right subunit", "Left subunit"))
legend("topright", legend=c("ATP-free", "ATP-bound"), cex=cex,
       fill=c("gray20", "gray90"))
#dev.off()

### plotting the different contributions. 
ylim=c(-4.5, 2.5)
s.inds <- ((diff.inds-1)%%524)+1
hmm <- NULL
hmm = t(cbind(apo$vdw,  holo$vdw, apo$ele, holo$ele,apo$sas, holo$sas))

#pdf("difference_residues_decomp.pdf", h=3, w=12)
#par(mfcol=c(1,1), mar=c(3,3,2,1), mgp=c(1.75,.6,0))
mp <- barplot(hmm[,diff.inds], beside=T, ylim=ylim, cex.lab=cex.label,
              ylab="Energy (kcal/mol)", cex.axis=cex.axis, 
              col=c("gray20","gray20","gray50","gray50", "gray90", "gray90"))
mtext(1, line=1, at=colMeans(mp), las=2, cex=cex,
      text=paste(seq[s.inds], s.inds+1))
mtext(3, at=0, text="B", cex=1.1)
legend("topright", legend=c("VDW", "ELE", "SAS"), cex=cex,
       fill=c("gray20", "gray50", "gray90"), ncol=3)
dev.off()











par(mfcol=c(2,1))
inds<-which(abs(a1$deltas2$total[,"total"])>0.5)

mp<-barplot(a1$deltas2$total[,"total"])
text(mp[inds], a1$deltas2$total[inds,"total"], labels=paste(s[inds], inds+1))

inds<-which(abs(b1$deltas2$total[,"total"])>0.5)
mp<-barplot(b1$deltas2$total[,"total"])
text(mp[inds], b1$deltas2$total[inds,"total"], labels=paste(s[inds], inds+1))




## which residues are the peaks
l.inds <- which(apo$tots[1:524]>1.5)
s.inds <- which(holo$tots[1:524]<(-1.5))


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


## print the whole shait
pdf("all_residues.pdf", w=12, h=8)
par(mfcol=c(4,1), mar=c(3,3,1,1), mgp=c(2,0.75,0))

inds <- c(1:524)
l.inds=c(33,50,74,215,389)
s.inds=c(38,46,182,269,385)

mp <- barplot(apo$total[inds], ylim=ylim, ylab=ylab)
mtext(3, at=-20, adj=0, text="A", cex=1.25, line=-1)
text(mp[s.inds], apo$total[s.inds]-offset, labels=paste(s[s.inds], s.inds+1, sep=""), las=2, cex=label.cex)
text(mp[l.inds], apo$total[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$total[inds], ylim=ylim, ylab=ylab)
mtext(3, at=-20, adj=0, text="B", cex=1.25, line=-1)
text(mp[s.inds], holo$total[s.inds]-offset, labels=paste(s[s.inds], s.inds+1, sep=""), las=2, cex=label.cex)
text(mp[l.inds], holo$total[l.inds]+offset, labels=paste(s[l.inds], l.inds+1, sep=""), las=2, cex=label.cex)


inds <- c(525:1048)
l.inds <- which(apo$total[inds]>1)
s.inds <- which(holo$total[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$total[inds], ylim=ylim, ylab=ylab)
mtext(3, at=-20, adj=0, text="C", cex=1.25, line=-1)
text(mp[s.inds], apo$total[s.inds+524]-offset, labels=paste(s[s.inds], s.inds+1, sep=""), las=2, cex=label.cex)
text(mp[l.inds], apo$total[l.inds+524]+offset, labels=paste(s[l.inds], l.inds+1, sep=""), las=2, cex=label.cex)

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



## zooming in on a area
## with errorbars now
n <- 2000*4
ylim=c(-6,3.5)
diff.inds <- which(abs(apo$total-holo$total)>0.5)

pdf("diff_resi.pdf", w=10, h=6)
mp <- plot.zoom(apo, holo, diff.inds, n, seq, ylim)
#mtext(1, line=3, at=c(mean(mp[,1:8]), mean(mp[,9:15])),
#      text=c("Right subunit", "Left subunit"))
legend("topright", legend=c("ATP-free", "ATP-bound"),
       fill=c("gray20", "gray90"))
dev.off()


"plot.zoom" <- function(apo, holo, inds, n, seq, ylim=NULL, ylab="Energy (kcal/mol)") {
  s.inds <- ((inds-1)%%524)+1
  hmm <- matrix(c(apo$total[inds], holo$total[inds]), nrow=2, byrow=T)
  print(head(hmm))
  mp <- barplot(hmm, ylim=ylim, ylab=ylab, beside=T)

  stderr=(apo$totals.stds[inds]/sqrt(n))
  errbar(mp[1,], apo$total[inds], apo$total[inds]+stderr, apo$total[inds]-stderr, add=T, cex=0.65)

  stderr=(holo$totals.stds[inds]/sqrt(n))
  errbar(mp[2,], holo$total[inds], holo$total[inds]+stderr, holo$total[inds]-stderr, add=T, cex=0.65)
  mtext(1, at=colMeans(mp), text=paste(seq[s.inds], s.inds+1), las=2)

  return(mp)
}




### plotting the different contributions. 
ylim=c(-4, 2.5)
s.inds <- ((diff.inds-1)%%524)+1
hmm <- NULL
hmm = t(cbind(apo$vdw, apo$ele, apo$sas))

pdf("difference_residues_decomp.pdf", h=8, w=10)
par(mfcol=c(2,1), mar=c(3,3,2,1))
mp <- barplot(hmm[,diff.inds], beside=T, ylim=ylim)
mtext(1, line=1, at=colMeans(mp), las=2,
      text=paste(seq[s.inds], s.inds+1))
mtext(3, at=0, text="A (ATP-free)", cex=1.2)
legend("topright", legend=c("VDW", "ELE", "SAS"), 
       fill=c("gray20", "gray50", "gray90"))


hmm <- NULL
hmm = t(cbind(holo$vdw, holo$ele, holo$sas))

mp <- barplot(hmm[,diff.inds], beside=T, ylim=ylim)
mtext(1, line=3, at=colMeans(mp), las=2, 
      text=paste(seq[s.inds], s.inds+1))
mtext(3, at=0, text="B (ATP-bound)", cex=1.2)
dev.off()




### which resiudes contribute most - APO only

ylim=c(-6, 2.5)
inds <- which(abs(apo$total)>1)
s.inds <- ((inds-1)%%524)+1
hmm <- NULL
hmm = t(cbind(apo$vdw, apo$ele, apo$sas))

pdf("residues_largest_contrib.pdf", h=6, w=10)
par(mfcol=c(1,1), mar=c(3,3,2,1))
mp <- barplot(hmm[,inds], beside=T, ylim=ylim)
mtext(1, line=1, at=colMeans(mp), las=2,
      text=paste(seq[s.inds], s.inds+1), cex=0.8)
#mtext(3, at=0, text="A (ATP-free)", cex=1.2)
legend("topright", legend=c("VDW", "ELE", "SAS"), 
       fill=c("gray20", "gray50", "gray90"))
dev.off()





















par(mfcol=c(2,1))
vdw <- as.numeric( a1$deltas2$total[diff.inds, "vdw"] )
ele <- as.numeric( a1$deltas2$total[diff.inds, c("ele", "polsol")] )
ele <- rowMeans(matrix(ele, ncol=2))
sas <- as.numeric( a1$deltas2$total[diff.inds, "nonpolsol"] )

hmm <- NULL
hmm = t(cbind(vdw, ele, sas))

mp <- barplot(hmm, beside=T)
mtext(1, line=3, at=colMeans(mp), las=2,
      text=paste(seq[s.inds], s.inds+1))


vdw <- as.numeric( b1$deltas2$total[diff.inds, "vdw"] )
ele <- as.numeric( b1$deltas2$total[diff.inds, c("ele", "polsol")] )
ele <- rowMeans(matrix(ele, ncol=2))
sas <- as.numeric( b1$deltas2$total[diff.inds, "nonpolsol"] )

hmm <- NULL
hmm = t(cbind(vdw, ele, sas))

mp <- barplot(hmm, beside=T)
mtext(1, line=3, at=colMeans(mp), las=2,
      text=paste(seq[s.inds], s.inds+1))


#legend("topright", legend=c("ATP-free", "ATP-bound"),
#       fill=c("gray20", "gray90"))








tot <- as.numeric( a1$deltas$total[diff.inds, "total"] )
std <- as.numeric( a1$deltas$total[diff.inds, 15] )
std = std / sqrt(500)
mp <- barplot(tot)
errbar(mp[,1], tot, std+tot, tot-std, add=T)




#mp <- barplot(a1$deltas2$total[, "total"], ylim=ylim, ylab=ylab)









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)