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

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

a <- read.mmpbsa.decomp("116/concat_FINAL_DECOMP_MMPBSA_1.dat")
b <- read.mmpbsa.decomp("129/concat_FINAL_DECOMP_MMPBSA_1.dat")

#a10ns <- read.mmpbsa.decomp("116/concat10n_FINAL_DECOMP_MMPBSA_1.dat")
#b10ns <- read.mmpbsa.decomp("129/concat10n_FINAL_DECOMP_MMPBSA_1.dat")

  
apo <- calc.means(list(a), n=1)
holo <- calc.means(list(b), n=1)

#apo <- calc.means(list(a10ns), n=1)
#holo <- calc.means(list(b10ns), n=1)


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

par(mfcol=c(2,1))
mp <- barplot(a$deltas2$total[inds,"total"], ylim=ylim, ylab=ylab)
text(mp[s.inds], a$deltas2$total[s.inds, "total"]-offset, labels=paste(s[s.inds], s.inds+1, sep=""), las=2, cex=label.cex)
text(mp[l.inds], a$deltas2$total[l.inds, "total"], labels=paste(s[l.inds], l.inds+1, sep=""), las=2, cex=label.cex)


mp <- barplot(b$deltas2$total[inds,"total"], ylim=ylim, ylab=ylab)
text(mp[s.inds], b$deltas2$total[s.inds, "total"]-offset, labels=paste(s[s.inds], s.inds+1, sep=""), las=2, cex=label.cex)
text(mp[l.inds], b$deltas2$total[l.inds, "total"], labels=paste(s[l.inds], l.inds+1, sep=""), las=2, cex=label.cex)

inds <- c(525:1048)
l.inds=c(12,116,224,251,302)
s.inds=c(7,68,113,196,230,280,284,359,518)

par(mfcol=c(2,1))
mp <- barplot(a$deltas2$total[inds,"total"], ylim=ylim, ylab=ylab)
text(mp[s.inds], a$deltas2$total[s.inds+524, "total"]-offset, labels=paste(s[s.inds], s.inds+1, sep=""), las=2, cex=label.cex)
text(mp[l.inds], a$deltas2$total[l.inds+524, "total"], labels=paste(s[l.inds], l.inds+1, sep=""), las=2, cex=label.cex)


mp <- barplot(b$deltas2$total[inds,"total"], ylim=ylim, ylab=ylab)
text(mp[s.inds], b$deltas2$total[s.inds+524, "total"]-offset, labels=paste(s[s.inds], s.inds+1, sep=""), las=2, cex=label.cex)
text(mp[l.inds], b$deltas2$total[l.inds+524, "total"], labels=paste(s[l.inds], l.inds+1, sep=""), las=2, cex=label.cex)
}



cex=1.4
cex.axis=1.4
cex.label=1.4
n <- 3600
ylim=c(-5.5,3.5)
diff.inds <- which(abs(apo$total-holo$total)>.50)
grps=c(10,9)
grps=cumsum(grps)


pdf("decomp_desiduess.pdf", w=8, h=10)
par(mfcol=c(2,1), mar=c(3,3,2,1), mgp=c(1.75,.5,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=2.5)

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 (right)")
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="cis (left)")
lines(x=c(mp[1,a], mp[2,b]), y=y)


#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. 
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)", las=2, cex.axis=cex.axis, 
              col=c("gray20","gray20","gray50","gray50", "gray90", "gray90"))
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 (right)")
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="cis (left)")
lines(x=c(mp[1,a], mp[nrow(mp),b]), y=y)

mtext(1, line=-.2, at=colMeans(mp), las=2, cex=cex,
      text=paste(seq[s.inds], s.inds+1))
mtext(3, at=0, text="B", cex=2.5)
legend("topright", legend=c("VDW", "ELE", "SAS"), cex=cex,
       fill=c("gray20", "gray50", "gray90"), ncol=3)
dev.off()



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

mp <- barplot(hmm[,diff.inds], beside=T, ylim=ylim,
              ylab="Energy (kcal/mol)")
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)
dev.off()





## write to file
inds <- 1:1048
s.inds <- ((c(1:1048)-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)


b <- NULL
for ( i in 1:length(pdb.dimer$atom[,"resno"])) {
  resno=as.numeric(pdb.dimer$atom[i,"resno"])
  b=c(b, apo$total[resno])
}
write.pdb(pdb.dimer, b=b)  







a1 <- read.mmpbsa.decomp("116/FINAL_DECOMP_MMPBSA_1.dat")
a2 <- read.mmpbsa.decomp("116/FINAL_DECOMP_MMPBSA_2.dat")
a3 <- read.mmpbsa.decomp("116/FINAL_DECOMP_MMPBSA_3.dat")
a4 <- read.mmpbsa.decomp("116/FINAL_DECOMP_MMPBSA_4.dat")
a5 <- read.mmpbsa.decomp("116/FINAL_DECOMP_MMPBSA_5.dat")
#a6 <- read.mmpbsa.decomp("116/FINAL_DECOMP_MMPBSA_6.dat")

b1 <- read.mmpbsa.decomp("129/FINAL_DECOMP_MMPBSA_1.dat")
b2 <- read.mmpbsa.decomp("129/FINAL_DECOMP_MMPBSA_2.dat")
b3 <- read.mmpbsa.decomp("129/FINAL_DECOMP_MMPBSA_3.dat")
b4 <- read.mmpbsa.decomp("129/FINAL_DECOMP_MMPBSA_4.dat")
b5 <- read.mmpbsa.decomp("129/FINAL_DECOMP_MMPBSA_5.dat")
#b6 <- read.mmpbsa.decomp("129/FINAL_DECOMP_MMPBSA_6.dat")


#par(mfcol=c(3,1))
#mp <- barplot(a1$deltas2$total[,"total"], ylim=ylim, ylab=ylab)
#mp <- barplot(a5$deltas2$total[,"total"], ylim=ylim, ylab=ylab)
#mp <- barplot(a6$deltas2$total[,"total"], ylim=ylim, ylab=ylab)


apo <- calc.means(list(a1,a2,a3,a4,a5), n=5)
holo <- calc.means(list(b1,b2,b3,b4,b5), n=5)

## 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(-9,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(2,3,0.5,1), mgp=c(2,0.75,0))

inds <- c(1:524)
l.inds=c(33,50,74,215,389)
s.inds=c(38,46,59,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, 482)
s.inds=c(38,46,60,182,241,269,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(3,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()








### which resiudes contribute most - APO only

ylim=c(-8, 2.5)
inds <- which(abs(apo$total)>1.5)
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"), cex=0.8,
       fill=c("gray20", "gray50", "gray90"))
dev.off()





## differences in the complex... (not deltaG)

pdf("internal_monomeric_diffs.pdf", w=15, h=8)
par(mfcol=c(1,1), mar=c(2,2,1,1))
inds <- c(1:524)
diff<-as.numeric(a10ns$ligand$total[inds,"total"])-as.numeric(b10ns$ligand$total[inds,"total"])
inds<-which(abs(diff)>1)
s.inds <- ((inds-1)%%524)+1
mp<-barplot(diff, ylim=c(-6,18))
text(mp[inds], diff[inds], labels=paste(s[s.inds], s.inds+1))
text(5,10, "favorable in ATP-bound", adj=0)
text(5,-4, "favorable in ATP-free", adj=0)
mtext(3, line=-1, at=0, text="Monomer A", cex=1.5)
dev.off()





diff.inds=inds


par(mfcol=c(2,1))
vdw1 <- as.numeric( a10ns$ligand$total[diff.inds, "vdw"] )
ele1 <- as.numeric( a10ns$ligand$total[diff.inds, c("ele", "polsol")] )
ele1 <- rowMeans(matrix(ele1, ncol=2))
sas1 <- as.numeric( a10ns$ligand$total[diff.inds, "nonpolsol"] )

hmm1 <- NULL
hmm1 = t(cbind(vdw1, ele1, sas1))

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


vdw2 <- as.numeric( b10ns$ligand$total[diff.inds, "vdw"] )
ele2 <- as.numeric( b10ns$ligand$total[diff.inds, c("ele", "polsol")] )
ele2 <- rowMeans(matrix(ele2, ncol=2))
sas2 <- as.numeric( b10ns$ligand$total[diff.inds, "nonpolsol"] )

hmm2 <- NULL
hmm2 = t(cbind(vdw2, ele2, sas2))

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



hmm <- hmm1-hmm2
mp <- barplot(hmm, beside=T)
mtext(1, line=0, 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)








}