library(bio3d)
source("mmpbsa_funs.R")

pdb <- read.pdb("1XCK_chainA_noWAT_noH.pdb")
#pdb.dimer <- read.pdb("complex.pdb")
seq <- seq.pdb(pdb)
s <- array(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")
c1 <- read.mmpbsa.decomp("129_noATP/FINAL_DECOMP_MMPBSA_1.dat")

                                     
#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



d<-a1$deltas$total-c1$deltas$total
d.inds <- which(abs(d)>1)

pdf("diff_mmpbsa_decomp.pdf", w=12, h=6)
par(mfcol=c(3,1), mar=c(2,2,1,1))
plot(a1$deltas$total, type="o", cex=0.5)
plot(c1$deltas$total, type="o", cex=0.5)
#plot(c1$deltas$total, type="o", cex=0.5)
plot(d, type="o", cex=0.5, col="red")
text(d.inds, d[d.inds]*1.05, labels=paste(s[d.inds], d.inds+1), col="red")
dev.off()



par(mfcol=c(2,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)