library(maptools)
library(rparallel)
library(bio3d)
library(ncdf)
library(plotrix)
source("dmat_funs.R")

##pdb <- read.pdb("1XCK_reference.pdb")
pdb <- read.pdb("1XCK_chainA_noWAT_noH.pdb")
pdb.ref <- pdb


seq <- seq.pdb(pdb)
s <- array(seq)
sse <- dssp(pdb)

ca.inds <- atom.select(pdb, "calpha")
dim <- length(ca.inds$atom)

## whole shit
trj.apo <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/121_1XCK_chainA_apo/results/traj/10-40ns_short_600frames_noWAT_noH.nc")
trj.holo <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/126_1XCK_chainA_MGATP/results/traj/10-40ns_short_600frames_noWAT_noH.nc")

dm.apo <- dm.xyz.trj( trj.apo, pdb,  threads=7 )
dm.holo <- dm.xyz.trj( trj.holo, pdb, threads=7 )



## Plots elements closer than 4Å in 50% of the frames
pdf("cmap.pdf")
g<-which(dm.apo>0.5, arr.ind=TRUE)
plot.bio3d(g[,1], g[,2], type='p', pch=2, col='green', sse=sse, top=FALSE, sse.border=TRUE)

g<-which(dm.holo>0.5, arr.ind=TRUE)
lines(g[,1], g[,2], type='p', pch=13, col='red')
dev.off()



c.md <- dcm(dm.apo, dm.holo, dim)

pdf("dcm_apo_holo_eq.pdf", height=12)
par(mfcol=c(2,1))
plot.dcm(c.md, pdb, xlim=c(1,135), ylim=c(1,135), legend=c('Apo', 'Holo'), legend.xy=c(100,40))
plot.dcm(c.md, pdb, xlim=c(405,525), ylim=c(405,525), legend=c('Apo', 'Holo'), legend.xy=c(500,440))
dev.off()


######################################

library(maptools)
library(rparallel)
library(bio3d)
library(ncdf)
library(plotrix)

source("dmat_funs.R")


pdb <- read.pdb("1XCK_chainA_noWAT_noH.pdb")


seq <- seq.pdb(pdb)
s <- array(seq)

ca.inds <- atom.select(pdb, "calpha")
dim <- length(ca.inds$atom)

#dir.prefix <- "/net/lutefisk/slars/groel_md/1XCK_chainA/"
#dir.prefix <- "~/lutefisk/groel_md/1XCK_chainA/"
#trj.suffix <- "/results/traj/10-40ns_noWAT_noH_noATP_noH.nc"

#apo.files <- c("121_1XCK_chainA_apo", "121_1XCK_chainA_apo_p1", "121_1XCK_chainA_apo_p2", "121_1XCK_chainA_apo_p3", "121_1XCK_chainA_apo_p4", "121_1XCK_chainA_apo_p5")
#holo.files <- c("126_1XCK_chainA_MGATP", "126_1XCK_MGATP_p1", "126_1XCK_MGATP_p2", "126_1XCK_MGATP_p3", "126_1XCK_MGATP_p4", "126_1XCK_MGATP_p5")

dir.prefix <- "/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/"
apo.files <- c("100ns_noWAT_noATP_noH_1000frames_p1.nc", "100ns_noWAT_noATP_noH_1000frames_p2.nc", "100ns_noWAT_noATP_noH_1000frames_p3.nc", "100ns_noWAT_noATP_noH_1000frames_p4.nc")
holo.files <- c("100ns_noWAT_noATP_noH_1000frames_p5.nc", "100ns_noWAT_noATP_noH_1000frames_p6.nc", "100ns_noWAT_noATP_noH_1000frames_p7.nc", "100ns_noWAT_noATP_noH_1000frames_p8.nc")


trj.apo <- NULL
#read.ncdf( paste(dir.prefix, apo.files[1], sep="") )
trj.holo <- NULL
#read.ncdf( paste(dir.prefix, holo.files[1], sep="") )

for ( i in 1:3) {
  trj <- read.ncdf(paste(dir.prefix, apo.files[i], sep=""))
  trj.apo <- rbind(trj.apo, trj[501:1000,])
}
trj <- read.ncdf(paste(dir.prefix, apo.files[4], sep=""))
trj.apo <- rbind(trj.apo, trj[501:973,])


for ( i in 1:length(holo.files)) {
  trj <- read.ncdf(paste(dir.prefix, holo.files[i], sep=""))
  trj.holo <- rbind(trj.holo, trj[501:1000,])
}


threads<-3
dm.apo <- dm.xyz.trj( trj.apo, pdb,  threads=threads )
save(dm.apo.50ns, file="dmat_apo_last50ns_2000frames.RData")

dm.holo <- dm.xyz.trj( trj.holo, pdb,  threads=threads )
save(dm.holo.50ns, file="dmat_holo_last50ns_2000frames.RData")


#save(dm.apo, dm.holo, file="dmat_apo_holo_600frames.RData")
#save(dm.apo, dm.holo, file="dmat_apo_holo_100ns_4000frames.RData")
##load("dmat_analysis.RData")
#load("dmat_apo_holo_3600frames.RData")

load("dmat_apo_100ns_4000frames.RData")
load("dmat_holo_100ns_4000frames.RData")

c.md <- dcm(dm.apo, dm.holo, occupancy=0.5, diff.cut=0.5)


pdf("dcm_apo_holo_MD.pdf", height=12, width=12)
par(mfcol=c(2,2))

source("dmat_funs.R")

xlim=c(1,140); ylim=c(405,524);
mat <- trunc.mat(c.md,xlim,ylim)
plot.dcm(mat, pdb,  xlim=xlim, ylim=ylim, legend=c('Apo', 'Holo'), legend.xy=c(100,40))

xlim=c(1,140); ylim=c(1,140);
mat <- trunc.mat(c.md,xlim,ylim)
plot.dcm(mat, pdb, xlim=xlim, ylim=ylim) 

xlim=c(405,524); ylim=c(405,524);
mat <- trunc.mat(c.md,xlim,ylim)
plot.dcm(mat, pdb, xlim=xlim, ylim=ylim) 

xlim=c(1,524); ylim=c(1,524);
mat <- trunc.mat(c.md,xlim,ylim)
plot.dcm(mat, pdb, legend=c('Apo, not holo', 'Holo, not apo'), legend.xy=c(400,200))


dev.off()







######################################


library(Hmisc)
library(maptools)
library(bio3d)
load("myresults2.RData")
source("dmat_funs.R")

#pdb.ref <- read.pdb("1XCK_chainA_noWAT_noH.pdb")
#pdb.ids <- pdbs$id[-exclude.inds]

#apo.ids <- pdb.ids[grep("1XCK_[ABCDEFGHIJKLMN]", pdb.ids)]
#holo.ids <- pdb.ids[grep("1SX4_[ABCDEFG]|1SVT_[ABCDEFG]", pdb.ids)]
#apo.ids <- pdb.ids[grep("1XCK_[A]", pdb.ids)]
#holo.ids <- pdb.ids[grep("1XCK_[ABC]", pdb.ids)]
#holo.ids <- pdb.ids[grep("1XCK_[I]", pdb.ids)]


#ids <- grep("1XCK_[ABC]|1KP8_[AB]", pdbs$id)
#apo.xray.dm <- dm.xyz.pdbs2(pdbs, ids, scut=4, dcut=4)


#ids <- grep("1SVT_[ABC]|1SX4_[AB]", pdbs$id)
#holo.xray.dm <- dm.xyz.pdbs2(pdbs, ids, scut=4, dcut=4)

#raw_pdbs_new/split_chain/3E76_F.pdb
#ids <- grep("1KP8_H|1XCK_A", pdbs$id)


#r.inds <- c(grep("1GRU_[HIJKLMN]", pdbs$id), grep("2C7E_[HIJKLMN]", pdbs$id))
#r1.inds <- c(grep("2C7E_[ABCDEFG]", pdbs$id))


## choosing only high-res x-ray structures

t.inds <- grep("1XCK|1SS8|1OEL", pdbs$id)
atp.inds <- grep("1KP8|1SX3", pdbs$id)
r2.inds <- grep("1AON_[ABCDEFG]|1SX4_[ABCDEFG]|1PF9_[ABCDEFG]|1SVT_[ABCDEFG]", pdbs$id)

library(rparallel)

t.xray.dm <- dm.xyz.pdbs(pdbs, t.inds, scut=4, dcut=4, threads=3)
atp.xray.dm <- dm.xyz.pdbs(pdbs, atp.inds, scut=4, dcut=4, threads=3)
r2.xray.dm <- dm.xyz.pdbs(pdbs, r2.inds, scut=4, dcut=4, threads=3)
#r1.xray.dm <- dm.xyz.pdbs(pdbs, r1.inds, scut=4, dcut=4, threads=2)
#r1.xray.dm <- dm.xyz.pdbs(pdbs, r1.inds, scut=4, dcut=4, threads=2)


dcm.tr2.xray <- dcm(t.xray.dm, r2.xray.dm, occupancy=0.5, diff.cut=0.5)
dcm.tatp.xray <- dcm(t.xray.dm, atp.xray.dm, occupancy=0.5, diff.cut=0.5)


#apo.xray.dm <- dm.xyz.pdbs(apo.ids, scut=4, dcut=4)
#holo.xray.dm <- dm.xyz.pdbs(holo.ids, scut=4, dcut=4)
##save(apo.xray.dm, holo.xray.dm, file="apo_holo_xray-DM.Rdata")
#save(t.xray.dm, r2.xray.dm, file="apo_holo_xray-DM_ALLstructures.Rdata")
#save(t.xray.dm, atp.xray.dm, r2.xray.dm, dcm.tr2.xray, dcm.tatp.xray, file="xray_DCM.RData"  )



xlim=c(191,376); ylim=c(191,376);
pdf("dcm_tmp.pdf")
plot.dcm(c.xray, pdb,  xlim=xlim, ylim=ylim)
dev.off()


pdf("dcm_apo_holo_xrays.pdf", height=12, width=12)
par(mfcol=c(2,2))

source("dmat_funs.R")

xlim=c(1,140); ylim=c(405,524);
mat <- trunc.mat(c.xray,xlim,ylim)
plot.dcm(mat, pdb.ref,  xlim=xlim, ylim=ylim, legend=c('Apo', 'Holo'), legend.xy=c(100,40))

xlim=c(1,140); ylim=c(1,140);
mat <- trunc.mat(c.xray,xlim,ylim)
plot.dcm(mat, pdb.ref, xlim=xlim, ylim=ylim) 

xlim=c(405,524); ylim=c(405,524);
mat <- trunc.mat(c.xray,xlim,ylim)
plot.dcm(mat, pdb.ref, xlim=xlim, ylim=ylim)

xlim=c(1,524); ylim=c(1,524);
plot.dcm(c.xray, pdb.ref,  xlim=xlim, ylim=ylim)


dev.off()




pdf("dcm_apo_holo_xrays_allInOne.pdf")
plot.dcm(c.xray, pdb.ref, legend=c('Apo, not holo', 'Holo, not apo'), legend.xy=c(400,200))
dev.off()





apo.xray.dm <- dm.xyz.pdbs("1XCK_chainA_noWAT_noH.pdb", scut=4, dcut=4)
apo.xray.dm2 <- dm.xyz.pdbs(apo.ids[1], scut=4, dcut=4)
c.xray <- dcm(apo.xray.dm, apo.xray.dm2, occupancy=0.5)
which(c.xray>0)




######################

library(plotrix)

source("dmat_funs.R")
#load("apo_holo_xray-DM_ALLstructures.Rdata")
load("xray_DCM.RData")
load("xrayPCA.RData")
load("dmat_apo_last50ns_2000frames.RData")
load("dmat_holo_last50ns_2000frames.RData")

pdb <- read.pdb("raw_pdbs/split_chain/1XCK_A.pdb")
pdb.ref <- pdb
seq <- seq.pdb(pdb)
s <- array(seq)
#sse <- stride(pdb)
sse <- dssp(pdb)


## contact difference contact maps
dcm.tr2.xray <- dcm(t.xray.dm, r2.xray.dm, occupancy=.5, diff.cut=.5)
dcm.tatp.xray <- dcm(t.xray.dm, atp.xray.dm, occupancy=.5, diff.cut=.5)
c.md <- dcm(dm.apo, dm.holo, occupancy=.5, diff.cut=.2)


## find common residues
sim <- which( (dcm.tr2.xray!=0 | dcm.tatp.xray!=0) & c.md!=0, arr.ind=T)
sim.r2 <- which(dcm.tr2.xray!=0 & c.md!=0, arr.ind=T)
sim.atp <- which(dcm.tatp.xray!=0 & c.md!=0, arr.ind=T)

## find difference in average distance
xray.avg.diff <- t.xray.dm$dmat.avg - r2.xray.dm$dmat.avg
xray.avg.diff2 <- t.xray.dm$dmat.avg - atp.xray.dm$dmat.avg
md.avg.diff <- dm.apo$dmat.avg - dm.holo$dmat.avg
##sim <- which(xray.avg.diff>dcut & md.avg.diff>dcut, arr.ind=T)
##sim <- which(dcm.tr2.xray!=0 & c.md!=0, arr.ind=T)


## list the residues found common in T/R2 and MD sim
sim.r2 <- cbind(sim.r2, round(xray.avg.diff[sim.r2],1), round(md.avg.diff[sim.r2],1))
diff <- data.frame(sim.r2)
diff = cbind(diff, (diff$V3+diff$V4)/2)
diff = cbind(diff, paste( s[diff[,"row"]], diff[,"row"]+1, sep="" ) )
diff = cbind(diff, paste( s[diff[,"col"]], diff[,"col"]+1, sep="" ) )

colnames(diff)<-c("row", "col", "avg.xray", "avg.md", "avg.diff", "res1", "res2")
diff.table1 <- diff[abs(order(diff$avg.diff)),]

neg.inds <- which(diff.table1[,"avg.xray"]<0 & diff.table1[,"avg.md"]<0, arr.ind=T)
pos.inds <- which(diff.table1[,"avg.xray"]>0 & diff.table1[,"avg.md"]>0, arr.ind=T)
paste( c(diff.table1[neg.inds,1]+1, diff.table1[neg.inds,2]+1), collapse="+", sep="" ) ## output for pymol
paste( c(diff.table1[pos.inds,1]+1, diff.table1[pos.inds,2]+1), collapse="+", sep="" ) ## output for pymol
diff.table1[which(diff.table1[,"avg.xray"]>0 & diff.table1[,"avg.md"]>0), c("res1", "res2", "avg.xray", "avg.md")]
diff.table1[which(diff.table1[,"avg.xray"]<0 & diff.table1[,"avg.md"]<0), c("res1", "res2", "avg.xray", "avg.md")]


## list the residues found common in T/R and MD sim
sim.atp <- cbind(sim.atp, round(xray.avg.diff[sim.atp],1), round(md.avg.diff[sim.atp],1))
diff <- data.frame(sim.atp)
diff = cbind(diff, (diff$V3+diff$V4)/2)
diff = cbind(diff, paste( s[diff[,"row"]], diff[,"row"]+1, sep="" ) )
diff = cbind(diff, paste( s[diff[,"col"]], diff[,"col"]+1, sep="" ) )

colnames(diff)<-c("row", "col", "avg.xray", "avg.md", "avg.diff", "res1", "res2")
diff.table2 <- diff[ order(abs(diff$avg.diff)),]
paste(c(diff.table2[, 1], diff.table2[, 2]), collapse="+") ## output for pymol

diff.table2[which(diff.table2[,"avg.xray"]>0 & diff.table2[,"avg.md"]>0), c("res1", "res2", "avg.xray", "avg.md")]
diff.table2[which(diff.table2[,"avg.xray"]<0 & diff.table2[,"avg.md"]<0), c("res1", "res2", "avg.xray", "avg.md")]


source("dcm_plot.r")






## these variabels is collected from the end of script in file dcm_plor.r
##sim
##res.labels.tr2
##res.labels.tatp
##int.conts$b


#tr2.inds <- which( c.md.xray!=0, arr.ind=T)
#which(tr2.inds==16)





md.inds <- which( c.md!=0, arr.ind=T)
fluctres <- c(16,27,28,33,51,52,53,55,58,60,62,86,420,474,477,479,480,484,135,136)


fluctres=paste("^", fluctres, "$", sep="")
fluctres.inds <- grep( paste(fluctres, collapse="|"), md.inds[,1])
md.avg <- round(md.avg.diff[md.inds[fluctres.inds,]], 1)
tr2.avg <- round( xray.avg.diff[md.inds[fluctres.inds,]],1)
tatp.avg <-  round(xray.avg.diff2[md.inds[fluctres.inds,]],1)
hoo <- cbind(md.inds[fluctres.inds,], md.avg, tr2.avg, tatp.avg)


fluctres.inds <- grep( paste(fluctres, collapse="|"), md.inds[,2])
md.avg.diff[md.inds[fluctres.inds,]]

md.avg <-  round(md.avg.diff[md.inds[fluctres.inds,]], 1)
tr2.avg <-  round(xray.avg.diff[md.inds[fluctres.inds,]], 1)
tatp.avg <-  round(xray.avg.diff2[md.inds[fluctres.inds,]], 1)
hoo <- rbind(hoo, cbind(md.inds[fluctres.inds,], md.avg, tr2.avg, tatp.avg))
hoo=data.frame(hoo)
hoo = hoo[ order(hoo$md.avg), ]


row.resi <- paste( s[hoo[,"row"]], hoo[,"row"]+1, sep="" )
col.resi <- paste( s[hoo[,"col"]], hoo[,"col"]+1, sep="" )
hoo=cbind(row.resi, col.resi, hoo)

res=paste("& ", hoo[, "row.resi"], hoo[,"col.resi"])
res=paste(res, hoo[, "tr2.avg"], hoo[,"md.avg"], hoo[,"tatp.avg"], sep="  &  ")
res=paste(res, "\\\\")

write(res, file="fluct_dcm.dat")
#cbind( res, hoo[, c("tr2.avg", "md.avg", "tatp.avg")] )
















for ( i in 1:nrow(sim) ) {
  a=sim[i,1]
  b=sim[i,2]

  tmp.inds <- which(int.conts$b[,"row"]==a | int.conts$b[,"row"]==b | 
                    int.conts$b[,"col"]==a | int.conts$b[,"col"]==b)

  print(int.conts$b[tmp.inds,])
}



sim.tmp <- int.conts$b[ which(int.conts$b[,"a.inds"]!="", arr.ind=T),]

a.inds=c()

for ( i in 1:nrow(sim.tmp))  {
  a=as.numeric(sim.tmp[i,"row"])
  b=as.numeric(sim.tmp[i,"col"])

  a.ind = as.numeric( unlist(strsplit(sim.tmp[i,"a.inds"], "|", fixed=T)) )
  a.inds=c(a.inds, a.ind)
}


int.conts$a[a.inds,]


for ( i in 1:nrow(sim) ) {
  a=sim[i,1]
  b=sim[i,2]

  tmp.inds <- which(int.conts$a[,"row"]==a | int.conts$a[,"row"]==b | 
                    int.conts$a[,"col"]==a | int.conts$a[,"col"]==b)

  print(paste(a, b))
  print(int.conts$a[tmp.inds,])
}




#diff = cbind(diff, (diff$V3+diff$V4)/2)
#colnames(diff)<-c("row", "col", "avg.xray", "avg.md", "avg.diff")
#diff[order(diff$avg.diff),]













#########







### ZZOOOM on Equatorial ###

pdf("dcm_apo_holo_xrays_md.pdf", height=12, width=12)
par(mfcol=c(3,3), mar=c(4,5,3,3))


xlim=c(1,100); ylim=c(20,120);
mat <- trunc.mat(dcm.tatp.xray,xlim,ylim)
plot.dcm(mat, pdb.ref, xlim=xlim, ylim=ylim,
         legend=c('Contacts Unbound','Contacts ATP Bound'), legend.xy=c(10,35),
         labels.cex=labels.cex, labels.offset=labels.offset,
         xlab="Residue No", ylab="Residue No", main="Xray Difference Contact Maps")
text(-10, 100, "A")
for ( i in 1:nrow(sim))  {
  draw.circle(sim[i,"row"], sim[i,"col"], 3, border="red", lty="dashed")
}


xlim=c(1,140); ylim=c(405,524);
mat <- trunc.mat(dcm.tatp.xray,xlim,ylim)
plot.dcm(mat, pdb.ref,  xlim=xlim, ylim=ylim,
         legend=c('Contacts Unbound','Contacts ATP Bound'), legend.xy=c(10,430),
         labels.cex=labels.cex, labels.offset=labels.offset,
         xlab="Residue No", ylab="Residue No")

for ( i in 1:nrow(sim))  {
  draw.circle(sim[i,"row"], sim[i,"col"], 3, border="red", lty="dashed")
}

xlim=c(405,500); ylim=c(440,510);
mat <- trunc.mat(dcm.tatp.xray,xlim,ylim)
plot.dcm(mat, pdb.ref, xlim=xlim, ylim=ylim,
         legend=c('Contacts Unbound','Contacts ATP Bound'), legend.xy=c(465,456),
         labels.cex=labels.cex, labels.offset=labels.offset,
         xlab="Residue No", ylab="Residue No")

for ( i in 1:nrow(sim))  {
  draw.circle(sim[i,"row"], sim[i,"col"], 3, border="red", lty="dashed")
}



xlim=c(1,100); ylim=c(20,120);
mat <- trunc.mat(dcm.tr2.xray,xlim,ylim)
plot.dcm(mat, pdb.ref, xlim=xlim, ylim=ylim,
         legend=c('Contacts Unbound','Contacts ATP Bound'), legend.xy=c(10,35),
         labels.cex=labels.cex, labels.offset=labels.offset,
         xlab="Residue No", ylab="Residue No", main="Xray Difference Contact Maps")
text(-10, 100, "A")
for ( i in 1:nrow(sim))  {
  draw.circle(sim[i,"row"], sim[i,"col"], 3, border="red", lty="dashed")
}


xlim=c(1,140); ylim=c(405,524);
mat <- trunc.mat(dcm.tr2.xray,xlim,ylim)
plot.dcm(mat, pdb.ref,  xlim=xlim, ylim=ylim,
         legend=c('Contacts Unbound','Contacts ATP Bound'), legend.xy=c(10,430),
         labels.cex=labels.cex, labels.offset=labels.offset,
         xlab="Residue No", ylab="Residue No")

for ( i in 1:nrow(sim))  {
  draw.circle(sim[i,"row"], sim[i,"col"], 3, border="red", lty="dashed")
}

xlim=c(405,500); ylim=c(440,510);
mat <- trunc.mat(dcm.tr2.xray,xlim,ylim)
plot.dcm(mat, pdb.ref, xlim=xlim, ylim=ylim,
         legend=c('Contacts Unbound','Contacts ATP Bound'), legend.xy=c(465,456),
         labels.cex=labels.cex, labels.offset=labels.offset,
         xlab="Residue No", ylab="Residue No")

for ( i in 1:nrow(sim))  {
  draw.circle(sim[i,"row"], sim[i,"col"], 3, border="red", lty="dashed")
}



xlim=c(1,100); ylim=c(20,120);
mat <- trunc.mat(c.md,xlim,ylim)
plot.dcm(mat, pdb.ref, xlim=xlim, ylim=ylim,
         legend=c('Contacts Unbound','Contacts ATP Bound'), legend.xy=c(10,35),
         labels.cex=labels.cex, labels.offset=labels.offset,
         xlab="Residue No", ylab="Residue No", main="MD  Difference Contact Maps")

for ( i in 1:nrow(sim))  {
  draw.circle(sim[i,"row"], sim[i,"col"], 3, border="red", lty="dashed")
}

xlim=c(1,140); ylim=c(405,524);
mat <- trunc.mat(c.md,xlim,ylim)
plot.dcm(mat, pdb.ref,  xlim=xlim, ylim=ylim,
         legend=c('Contacts Unbound','Contacts ATP Bound'), legend.xy=c(10,430),
         labels.cex=labels.cex, labels.offset=labels.offset,
         xlab="Residue No", ylab="Residue No")

for ( i in 1:nrow(sim))  {
  draw.circle(sim[i,"row"], sim[i,"col"], 3, border="red", lty="dashed")
}


xlim=c(405,500); ylim=c(440,510);
mat <- trunc.mat(c.md,xlim,ylim)
plot.dcm(mat, pdb.ref, xlim=xlim, ylim=ylim,
         legend=c('Contacts Unbound','Contacts ATP Bound'), legend.xy=c(465,456),
         labels.cex=labels.cex, labels.offset=labels.offset,
         xlab="Residue No", ylab="Residue No")

for ( i in 1:nrow(sim))  {
  draw.circle(sim[i,"row"], sim[i,"col"], 3, border="red", lty="dashed")
}

dev.off()





#############



source("dmat_funs.R")

c.xray <- dcm(apo.xray.dm, holo.xray.dm, occupancy=0.5, diff.cut=0.5)
c.md <- dcm(dm.apo, dm.holo, occupancy=0.5, diff.cut=0)

labels.cex=0.6
labels.offset=0




pdf("dcm_apo_holo_xrays_md_apical_intermediate.pdf", height=6, width=12)
par(mfcol=c(1,2))


xlim=c(130,410); ylim=c(130,410);
mat <- trunc.mat(c.xray,xlim,ylim)
plot.dcm(mat, pdb.ref, xlim=xlim, ylim=ylim,
         legend=c('Xray Apo-vs-Holo','Xray Holo-vs-Apo'), legend.xy=c(10,35),
         labels.cex=labels.cex, labels.offset=labels.offset )


xlim=c(130,410); ylim=c(130,410);
mat <- trunc.mat(c.md,xlim,ylim)
plot.dcm(mat, pdb.ref, xlim=xlim, ylim=ylim,
         legend=c('MD Apo-vs-Holo','MD Holo-vs-Apo'), legend.xy=c(10,35),
         labels.cex=labels.cex, labels.offset=labels.offset )

dev.off()










########## TRASH




## find important x-ray interactions
## ... those who change interaction partner...?
diff.inds <- which(dcm.tr2.xray!=0, arr.ind=T)
avg.dist.t <- t.xray.dm$dmat.avg[diff.inds] 
avg.dist.r2 <- r2.xray.dm$dmat.avg[diff.inds]
tr2.table = cbind(diff.inds, avg.dist.t, avg.dist.r2, avg.dist.t-avg.dist.r2)
colnames(tr2.table)<-c("row", "col", "avg.dist.t", "avg.dist.r2", "diff")
tr2.table=tr2.table[ order(tr2.table[,"diff"]), ]


t.inds <- which(tr2.table[,"diff"]<0)
r2.inds <- which(tr2.table[,"diff"]>0)

t <- tr2.table[t.inds,]
r2 <- tr2.table[r2.inds,]
flip <- c(1:nrow(r2))*0

#find change in interaction partners
r2.tmp <- r2
row.t <- c(1:nrow(r2.tmp))*0
col.t <- c(1:nrow(r2.tmp))*0
r2.tmp=cbind(r2.tmp, row.t, col.t)
r2=cbind(r2, flip)
tr2.changed.inds <- c()


for ( i in 1:nrow(r2.tmp) ) {

  resa <- r2[i,"row"]
  resb <- r2[i,"col"]

  row.t=c()
  col.t=c()

  j=which( t[,"row"]==resa | t[,"row"]==resb |
    t[,"col"]==resa | t[,"col"]==resb, arr.ind=T )
  tr2.changed.inds=c(tr2.changed.inds, j)
  
  j=which(t[,"row"]==resa)
  if(length(j)>0)
    col.t=c(col.t, t[j,"col"])

  j=which(t[,"col"]==resa)
  if(length(j)>0)
    row.t=c(row.t, t[j,"row"])

  j=which(t[,"row"]==resb)
  if(length(j)>0)
    col.t=c(col.t, t[j,"col"])
  
  j=which(t[,"col"]==resb)
  if(length(j)>0)
    row.t=c(row.t, t[j,"row"])

  r2.tmp[i, "row.t"]= paste(row.t, collapse="|")
  r2.tmp[i, "col.t"]= paste(col.t, collapse="|")
 
}

flip.inds <- c()
table=c()
for ( i in 1:nrow(r2.tmp) ) {

  if ( r2.tmp[i,"row.t"] != "" | r2.tmp[i,"col.t"] != "" ) {
    flip.inds=c(flip.inds,i)
    
    resa=r2.tmp[i,"row"]
    resb=r2.tmp[i,"col"]
    
    resa.t=r2.tmp[i,"row.t"]
    resb.t=r2.tmp[i,"col.t"]

    text=paste("R2-int: ", resa, resb, "... T-int: ",
      resa, resb.t, " || ", resb, resa.t)
    table=c(table, text)

  }
}
flip.inds <- which(r2.tmp[, "row.t"]!="" | r2.tmp[, "col.t"]!="" )
r2[flip.inds,"flip"]=1




## find important x-ray interactions
## ... those who change interaction partner...?
diff.inds <- which(dcm.tatp.xray!=0, arr.ind=T)
avg.dist.t <- t.xray.dm$dmat.avg[diff.inds] 
avg.dist.atp <- atp.xray.dm$dmat.avg[diff.inds]
tr.table = cbind(diff.inds, avg.dist.t, avg.dist.atp, avg.dist.t-avg.dist.atp)
colnames(tr.table)<-c("row", "col", "avg.dist.t", "avg.dist.atp", "diff")
tr.table=tr.table[ order(tr.table[,"diff"]), ]


t.inds <- which(tr.table[,"diff"]<0)
r.inds <- which(tr.table[,"diff"]>0)

t <- tr.table[t.inds,]
r <- tr.table[r.inds,]

#find change in interaction partners
row.t <- c(1:nrow(r))*0
col.t <- c(1:nrow(r))*0
t.inds <- c(1:nrow(r))*0
r.tmp <- r
r.tmp=cbind(r.tmp, row.t, col.t, t.inds)
r.tmp=data.frame(r.tmp)
r.tmp=transform(r.tmp, row.t=as.character(row.t),
                       col.t=as.character(col.t),
                       t.inds=as.character(t.inds) )

flip <- c(1:nrow(r))*0
r=cbind(r, flip)
tr.changed.inds <- c()

for ( i in 1:nrow(r.tmp) ) {

  resa <- r.tmp[i,"row"]
  resb <- r.tmp[i,"col"]

  row.t=c()
  col.t=c()

  k=which( t[,"row"]==resa | t[,"row"]==resb |
    t[,"col"]==resa | t[,"col"]==resb, arr.ind=T )
  tr.changed.inds=c(tr.changed.inds, k)

  j=which(t[,"row"]==resa)
  if(length(j)>0)
    col.t=c(col.t, t[j,"col"])

  j=which(t[,"col"]==resa)
  if(length(j)>0)
    row.t=c(row.t, t[j,"row"])

  j=which(t[,"row"]==resb)
  if(length(j)>0)
    col.t=c(col.t, t[j,"col"])
  
  j=which(t[,"col"]==resb)
  if(length(j)>0)
    row.t=c(row.t, t[j,"row"])

  r.tmp[i, "row.t"]= paste(row.t, collapse="|")
  r.tmp[i, "col.t"]= paste(col.t, collapse="|")
  r.tmp[i, "t.inds"]= paste(k, collapse="|")
  
}

flip.inds <- c()
table=c()
for ( i in 1:nrow(r.tmp) ) {

  if ( r.tmp[i,"row.t"] != "" | r.tmp[i,"col.t"] != "" ) {
    flip.inds=c(flip.inds,i)
    
    resa=r.tmp[i,"row"]
    resb=r.tmp[i,"col"]
    
    resa.t=r.tmp[i,"row.t"]
    resb.t=r.tmp[i,"col.t"]

    text=paste("R-int: ", resa, resb, "... T-int: ",
      resa, resb.t, " || ", resb, resa.t)
    table=c(table, text)
  }
}
flip.inds <- which(r.tmp[, "row.t"]!="" | r.tmp[, "col.t"]!="" )
r[flip.inds,"flip"]=1


#sim.tmp <- r.tmp[which(r.tmp[,"t.inds"]!="", arr.ind=T),]