dirname(system("pwd", intern = TRUE))
library(bio3d)

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


seq <- seq.pdb(pdb)
blast <- blast.pdb(seq)

##print(blast)
#pdf("blastplot.pdf", onefile=TRUE)
hits <- plot.blast(blast, cutoff=500 ) 
#dev.off()
head(hits)


for ( i in 1:length(hits$pdb.id) ) {
	id <- hits$hits[i, "pdb.id"]
	if ( substr(id,1,4) == "1WF4" ) {
          
          hits$hits[i, "pdb.id"][1] <-  substr(id,1,6)
	}
        
}


for ( i in 1:length(hits$pdb.id) ) {
	id <- hits$pdb.id[i]
	if ( substr(id,1,4) == "1WF4" ) {
          hits$pdb.id[i][1] <-  substr(id,1,6)
	}
        
}



unq.ids <- unique(substr(hits$pdb.id, 1, 4))




chains <- unique(pdb$atom[,"chain"])
#print(chains)


raw.files <- get.pdb(unq.ids, path = "raw_pdbs_new")


split.pdb(raw.files, path = "raw_pdbs_new/split_chain/")
## 1WF4 remains a problem due to small letters in the written PDB-filename
files <- paste("raw_pdbs_new/split_chain/", hits$pdb.id, ".pdb", sep = "")

## align the sequences
aln <- pdbaln(files)

##save(files, aln, pdbs,pdb, blast, file="xrayPCA.RData")

## Core superposition ##
pdbs <- read.fasta.pdb(aln, "", "", het2atom = T)
gaps.res <- gap.inspect(pdbs$ali)
exclude.inds <- which(gaps.res$row > 7)

## remove more suckers - these have 2 additional residues...
exclude.inds = c(exclude.inds, grep("1WF4", hits$pdb.id), grep("1WE3", hits$pdb.id), grep("3C9V", hits$pdb.id)) 

files.excl <- files[-exclude.inds]

## for obtaining specific pdb-codes
files.excl = grep("1AON_[HIJKLMN]|1SX4_[HIJKLMN]|1PF9_[HIJKLMN]|1SVT_[HIJKLMN]|1XCK|1KP8|2C7E", hits$pdb.id)
aln <- pdbaln(files.excl)
pdbs <- read.fasta.pdb(aln, "", "", het2atom = T)

pdb.codes <- unique(substr(pdbs$id, 26,29))
table(substr(pdbs$id, 26,29))


core <- core.find(pdbs)
xyz <- fit.xyz(fixed = pdbs$xyz[1, ], mobile = pdbs, fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", outpath = "core_fitlsq_new/", full.pdbs = TRUE, het2atom = TRUE)

gaps.pos <- gap.inspect(pdbs$xyz)
gaps.res <- gap.inspect(pdbs$ali)


## lets skip the N and C terminal residues - earlier 9-515 has been used
#pdb1 <- read.pdb("1XCK_reference.pdb")
#ca.inds.ref <- atom.select(pdb1, "///9:515///CA/")

#pc.xray <- pca.xyz(xyz[27:1545,gaps.pos$f.inds])
pc.xray <- pca.xyz(xyz[,gaps.pos$f.inds])
#inds.pca <- 27:1544
#pc.xray <- pca.xyz(xyz[,inds.pca])

#save(files, aln, pdbs, xyz, pdb, blast, hits, core, gaps.pos, gaps.res, pc.xray, file="xrayPCA2.RData")
#load("xrayPCA.RData")

plot(pc.xray, pch = 1 )


library(gclus)
library(maptools)
library(bio3d)

rd <- rmsd(xyz[, gaps.pos$f.inds])
dis <- as.dist(rd)
hc <- hclust(dis)
hc1 <- reorder.hclust(hc, dis)
id <- substr(basename(pdbs$id), 1, 6)
grps <- cutree(hc1, k=3)


highlight.ind <-  c(grep("1GRU_H", id), grep("1XCK_A", id),
                    grep("1SVT_A", id),
                    grep("2C7E_A", id), grep("2C7E_H", id),
                    grep("1GR5_A", id), grep("1KP8_A", id),
                    grep("3CAU_L", id), grep("2C7C_B", id) )

pca.scree <- plot.pca.scree(pc.xray, xlim=c(1,10))


pdf("xray_PCA.pdf")
plot(pc.xray, cex=0.5)
dev.off()



grps=1
pdf("xray_conformerplot_labelALL.pdf")
plot(pc.xray$z[,1], pc.xray$z[,2], col=grps+1, xlab = "PC1", ylab = "PC2", cex=0.8, pch=20)
text(pc.xray$z[, 1], pc.xray$z[, 2], id, cex=0.8, pos=c(4,1,2,1,4,4,4,1), col="grey50")

#pointLabel(pc.xray$z[which(grps==grp), 1], pc.xray$z[which(grps==grp), 2], id[which(grps==grp)], offset=40, cex=0.6)
abline(h=0, lty=2, col="grey60")
abline(v=0, lty=2, col="grey60")
dev.off()


library(gridBase)
library(plotrix)
library(car)

pdf("xray_conformerplot_2.pdf", h=6, w=6)
plot(pc.xray$z[,1], pc.xray$z[,2], col=grps+1, xlab = "PC1", ylab = "PC2", cex=0.8, pch=20, ylim=c(-45, 180))
#pointLabel(pc.xray$z[highlight.ind, 1], pc.xray$z[highlight.ind, 2], id[highlight.ind], offset=5, cex=0.8)
text(pc.xray$z[highlight.ind, 1], pc.xray$z[highlight.ind, 2], id[highlight.ind], cex=0.8, pos=c(4,1,2,1,4,4,4,1), col="grey50")
abline(h=0, lty=2, col="grey60")
abline(v=0, lty=2, col="grey60")

## cis R conf
r.inds <- c(grep("2C7E_[ABCDEFG]", id))
v1 <- as.vector(pc.xray$z[r.inds,1])
v2 <- as.vector(pc.xray$z[r.inds,2])
m <- matrix(c(v1,v2), ncol=2)
ellipse(c(mean(v1),mean(v2)), cov(m), 4, lty="dashed", col="grey80", center.pch=0, lwd=0.7)
text(mean(v1)-12, mean(v2)+8, "Cis r (D398A)", cex=0.9, col="black")


## trans R conf
r.inds <- c(grep("1GRU_[HIJKLMN]", id), grep("2C7E_[HIJKLMN]", id))
v1 <- as.vector(pc.xray$z[r.inds,1])
v2 <- as.vector(pc.xray$z[r.inds,2])
m <- matrix(c(v1,v2), ncol=2)
ellipse(c(mean(v1),mean(v2)), cov(m), 2, lty="dashed", col="grey80", center.pch=0, lwd=0.7)
text(mean(v1)+1, mean(v2)+5, "Trans r", cex=0.9, col="black")

## T conf
r2.inds <- which(grps==2)
r2.inds=r2.inds[-r.inds]
v1 <- as.vector(pc.xray$z[r2.inds,1])
v2 <- as.vector(pc.xray$z[r2.inds,2])
m <- matrix(c(v1,v2), ncol=2)
ellipse(c(mean(v1),mean(v2)), cov(m), 2, lty="dashed", col="grey80", center.pch=0, lwd=0.7)
text(mean(v1)-18, mean(v2)+28, "Cis/Trans t", cex=0.9, col="black")


## R'' conf
v1 <- as.vector(pc.xray$z[which(grps==1),1])
v2 <- as.vector(pc.xray$z[which(grps==1),2])
m <- matrix(c(v1,v2), ncol=2)
ellipse(c(mean(v1),mean(v2)), cov(m), 5, lty="dashed", col="grey80", center.pch=0, lwd=0.7)
text(mean(v1)-7, mean(v2)+25, "Cis r''", cex=0.9, col="black")



vp <- baseViewports()
pushViewport(vp$inner, vp$figure, vp$plot)
pushViewport(viewport(x=.55,y=0.55, width=.45, height=.45, just=c("left", "bottom")))
#par(fig=gridFIG(), new=T)
par(plt=gridPLT(), new=T, mgp=c(0,0.35,0))
#grid.rect(gp=gpar(lwd=3, col="red"))
#plot.pca.scree(pc.xray, xlim=c(1,10), ylab="", cex=0.2)

plot(pca.scree$percent[1:4], type="o", cex=0.4, xlim=c(1,4), xlab="", ylab="", axes=F, frame.plot=T)
text(c(1:3), pca.scree$percent[1:3], labels=round(pca.scree$cumv[1:3],1), pos=c(4,2,3), cex=.85, offset=.4)

axis(2, at=c(0,25,50,75,100), label=c(0,25,50,75,100), las=2, tck=-0.02, cex.axis=0.85 )
mtext("Proportion of Variance (%)", side=2, line=1.1, cex=0.9)

axis(1, at=c(1,2,3,4), label=c(1,2,3,4), tck=-0.02, cex.axis=0.75)
mtext("Eigenvalue Rank", side=1, line=1.1, cex=0.9)

popViewport(4)

dev.off()


#pca.scree$percent


#pdf("xray_scree.pdf")
pca.scree <- plot.pca.scree(pc.xray, xlim=c(1,10))
#dev.off()





pdb.closed <- read.pdb("1XCK_reference.pdb")
ca.inds <- atom.select(pdb.closed, "calpha")

trj121 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/121_1XCK_chainA_apo/results/traj/300ns_short_15000frames_noWAT.nc")
trj126 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/126_1XCK_chainA_MGATP/results/traj/300ns_short_15000frames_noWAT_noATP.nc")
trj3 <- read.ncdf("/net/gulrotkake/slars/groel_md/1SVT/123_1SVT_chainA_apo/results/traj/300ns_short_15000frames_noWAT.nc")
trj4 <- read.ncdf("/net/gulrotkake/slars/groel_md/1SVT/124_1SVT_chain_MGADP/results/traj/300ns_short_15000frames_noWAT_noATP.nc")


xyz_traj1 <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  trj121[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", outpath = "core_fitlsq/", full.pdbs = FALSE, het2atom = TRUE)
xyz_traj2 <- fit.xyz(fixed =  pdbs$xyz[1, ], mobile =  trj126[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, db.path = "", pdbext = "", outpath = "core_fitlsq/", full.pdbs = FALSE, het2atom = TRUE)
xyz_traj3 <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  trj3[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", outpath = "core_fitlsq/", full.pdbs = FALSE, het2atom = TRUE)
xyz_traj4 <- fit.xyz(fixed =  pdbs$xyz[1, ], mobile =  trj4[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, db.path = "", pdbext = "", outpath = "core_fitlsq/", full.pdbs = FALSE, het2atom = TRUE)


highlight.ind2 <-  c(grep("1XCK_A", id), grep("1SVT_A", id), grep("2C7E_A", id))

d121 <- pca.project(xyz_traj1[,gaps.pos$f.inds], pc.xray)
d126 <- pca.project(xyz_traj2[,gaps.pos$f.inds], pc.xray)
d123 <- pca.project(xyz_traj3[,gaps.pos$f.inds], pc.xray)
d124 <- pca.project(xyz_traj4[,gaps.pos$f.inds], pc.xray)



xlims=c(-150,300)
ylims=c(-50,250)

pdf("conformerplot_xray_md.pdf")
par(mfrow=c(2,2), mar=c(4,4,3,2))

plot( d121[,1],d121[,2], type="p", main = "Closed Unbound", xlab="PC1", ylab="PC2", col=densCols(d121[,1], d121[,2]), pch=20, xlim = xlims, ylim=ylims)
points( pc.xray$z[,1], pc.xray$z[,2], col = "black", pch = 16, cex = 1 )
points( pc.xray$z[,1], pc.xray$z[,2], col="gray", pch = 16, cex = 0.4 )
points(pc.xray$z[highlight.ind2, 1], pc.xray$z[highlight.ind2, 2], col = c("orange", "green", "red"), pch = 16 )
#text(-300+10, 250-10, "A", cex=1.5)
mtext("A", side=3, at=xlims[1], cex=1.5, line=1)


pdf("confplot_126.pdf", w=8, h=8)
plot( d126[,1],d126[,2], type="p", main = "", xlab="PC1", ylab="PC2", col=densCols(d126[,1], d126[,2]), pch=20, xlim = xlims, ylim=ylims)
points( pc.xray$z[,1], pc.xray$z[,2], col = "black", pch = 16, cex = 1 )
points( pc.xray$z[,1], pc.xray$z[,2], col="gray", pch = 16, cex = 0.4 )
points(pc.xray$z[highlight.ind2, 1], pc.xray$z[highlight.ind2, 2], col = c("orange", "green", "red"), pch = 16 )
#text(-300+10, 250-10, "B", cex=1.5)
#mtext("B", side=3, at=xlims[1], cex=1.5, line=1)
dev.off()


xlims=c(-150,370)
ylims=c(-120,250)

plot( d123[,1],d123[,2], type="p", main = "Open Unbound", xlab="PC1", ylab="PC2", col=densCols(d123[,1], d123[,2]), pch=20, xlim = xlims, ylim=ylims)
points( pc.xray$z[,1], pc.xray$z[,2], col = "black", pch = 16, cex = 1 )
points( pc.xray$z[,1], pc.xray$z[,2], col="gray", pch = 16, cex = 0.4 )
points(pc.xray$z[highlight.ind2, 1], pc.xray$z[highlight.ind2, 2], col = c("orange", "green", "red"), pch = 16 )
#text(-400+15, 175-10, "C", cex=1.5)
mtext("C", side=3, at=xlims[1], cex=1.5, line=1)

plot( d124[,1],d124[,2], type="p", main = "Open ADP Bound", xlab="PC1", ylab="PC2", col=densCols(d124[,1], d124[,2]), pch=20, xlim = xlims, ylim=ylims)
points( pc.xray$z[,1], pc.xray$z[,2], col = "black", pch = 16, cex = 1 )
points( pc.xray$z[,1], pc.xray$z[,2], col="gray", pch = 16, cex = 0.4 )
points(pc.xray$z[highlight.ind2, 1], pc.xray$z[highlight.ind2, 2], col = c("orange", "green", "red"), pch = 16 )
#text(-400+15, 175-10, "D", cex=1.5)
mtext("D", side=3, at=xlims[1], cex=1.5, line=1)

dev.off()


##pc.xray <- pca.xyz(xyz[,gaps.pos$f.inds])
pc1 <- mktrj.pca(pc.xray, pc=1, file="pc1.pdb")
pc2 <- mktrj.pca(pc.xray, pc=2, file="pc2.pdb")
pc3 <- mktrj.pca(pc.xray, pc=3, file="pc3.pdb")









trj1 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/100ns_noWAT_noATP_5000frames_p1.nc")
trj2 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/100ns_noWAT_noATP_5000frames_p2.nc")
trj3 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/100ns_noWAT_noATP_5000frames_p3.nc")
trj4 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/100ns_noWAT_noATP_5000frames_p4.nc")
trj5 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/100ns_noWAT_noATP_5000frames_p5.nc")
trj6 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/100ns_noWAT_noATP_5000frames_p6.nc")

trj7 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/100ns_noWAT_noATP_5000frames_p7.nc")
trj8 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/100ns_noWAT_noATP_5000frames_p8.nc")
trj9 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/100ns_noWAT_noATP_5000frames_p9.nc")
trj10 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/100ns_noWAT_noATP_5000frames_p10.nc")
trj11 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/100ns_noWAT_noATP_5000frames_p11.nc")
trj12 <- read.ncdf("/net/lutefisk/slars/groel_md/1XCK_chainA/auto_calcs/trajs/100ns_noWAT_noATP_5000frames_p12.nc")


xyz1 <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  trj1[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz2 <- fit.xyz(fixed =  pdbs$xyz[1, ], mobile =  trj2[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, db.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz3 <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  trj3[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz4 <- fit.xyz(fixed =  pdbs$xyz[1, ], mobile =  trj4[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, db.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz5 <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  trj5[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz6 <- fit.xyz(fixed =  pdbs$xyz[1, ], mobile =  trj6[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, db.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz7 <- fit.xyz(fixed = pdbs$xyz[1, ], mobile =  trj7[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, pdb.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz8 <- fit.xyz(fixed =  pdbs$xyz[1, ], mobile =  trj8[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, db.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz9 <- fit.xyz(fixed =  pdbs$xyz[1, ], mobile =  trj9[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, db.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz10 <- fit.xyz(fixed =  pdbs$xyz[1, ], mobile =  trj10[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, db.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz11 <- fit.xyz(fixed =  pdbs$xyz[1, ], mobile =  trj11[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, db.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)
xyz12 <- fit.xyz(fixed =  pdbs$xyz[1, ], mobile =  trj12[, ca.inds$xyz], fixed.inds = core$c1A.xyz, mobile.inds = core$c1A.xyz, db.path = "", pdbext = "", full.pdbs = FALSE, het2atom = TRUE)



highlight.ind2 <-  c(grep("1XCK_A", id), grep("1SVT_A", id), grep("2C7E_A", id))

p1a <- pca.project(xyz1[1:2500, gaps.pos$f.inds], pc.xray)
p2a <- pca.project(xyz2[1:2500, gaps.pos$f.inds], pc.xray)
p3a <- pca.project(xyz3[1:2500, gaps.pos$f.inds], pc.xray)
p4a <- pca.project(xyz4[1:2500, gaps.pos$f.inds], pc.xray)
p5a <- pca.project(xyz5[1:2500, gaps.pos$f.inds], pc.xray)
p6a <- pca.project(xyz6[1:2500, gaps.pos$f.inds], pc.xray)
p7a <- pca.project(xyz7[1:2500, gaps.pos$f.inds], pc.xray)
p8a <- pca.project(xyz8[1:2500, gaps.pos$f.inds], pc.xray)
p9a <- pca.project(xyz9[1:2500,gaps.pos$f.inds], pc.xray)
p10a <- pca.project(xyz10[1:2500,gaps.pos$f.inds], pc.xray)
p11a <- pca.project(xyz11[1:2500,gaps.pos$f.inds], pc.xray)
p12a <- pca.project(xyz12[1:2500,gaps.pos$f.inds], pc.xray)




p1 <- pca.project(xyz1[2501:nrow(xyz1), gaps.pos$f.inds], pc.xray)
p2 <- pca.project(xyz2[2501:nrow(xyz2), gaps.pos$f.inds], pc.xray)
p3 <- pca.project(xyz3[2501:nrow(xyz3), gaps.pos$f.inds], pc.xray)
p4 <- pca.project(xyz4[2501:nrow(xyz4), gaps.pos$f.inds], pc.xray)
p5 <- pca.project(xyz5[2501:nrow(xyz5), gaps.pos$f.inds], pc.xray)
p6 <- pca.project(xyz6[2501:nrow(xyz6), gaps.pos$f.inds], pc.xray)
p7 <- pca.project(xyz7[2501:nrow(xyz7), gaps.pos$f.inds], pc.xray)
p8 <- pca.project(xyz8[2501:nrow(xyz8), gaps.pos$f.inds], pc.xray)
p9 <- pca.project(xyz9[2501:nrow(xyz9),gaps.pos$f.inds], pc.xray)
p10 <- pca.project(xyz10[2501:nrow(xyz10),gaps.pos$f.inds], pc.xray)
p11 <- pca.project(xyz11[2501:nrow(xyz11),gaps.pos$f.inds], pc.xray)
p12 <- pca.project(xyz12[2501:nrow(xyz12),gaps.pos$f.inds], pc.xray)


p.apo <- round( c(p1[nrow(xyz1)-2501,1], p2[nrow(xyz2)-2501,1], p3[nrow(xyz3)-2501,1],
                  p4[nrow(xyz4)-2501,1], p5[nrow(xyz5)-2501,1], p6[nrow(xyz6)-2501,1] ), 1)
p.holo <- round( c(p7[nrow(xyz7)-2501,1], p8[nrow(xyz8)-2501,1], p9[nrow(xyz9)-2501,1],
                   p10[nrow(xyz10)-2501,1], p11[nrow(xyz11)-2501,1], p12[nrow(xyz12)-2501,1] ), 1)



m1<- round( c(mean(p1[,1]), mean(p2[,1]), mean(p3[,1]),
              mean(p4[,1]), mean(p5[,1]), mean(p6[,1])), 1)

m2<- round( c(mean(p7[,1]), mean(p8[,1]), mean(p9[,1]),
              mean(p10[,1]), mean(p11[,1]), mean(p12[,1])), 1)

med1<-round( c(median(p1[,1]), median(p2[,1]), median(p3[,1]),
               median(p4[,1]), median(p5[,1]), median(p6[,1]) ), 1)

med2<-round( c(median(p7[,1]), median(p8[,1]), median(p9[,1]),
               median(p10[,1]), median(p11[,1]), median(p12[,1]) ), 1)

var1<-round( c(var(p1[,1]), var(p2[,1]), var(p3[,1]),
               var(p4[,1]), var(p5[,1]), var(p6[,1]) ), 1)

var2<-round( c(var(p7[,1]), var(p8[,1]), var(p9[,1]),
               var(p10[,1]), var(p11[,1]), var(p12[,1]) ), 1)


sim.stats <- data.frame( m1, m2, var1, var2, med1, med2, p.apo, p.holo )
colnames(sim.stats) <- c("Mean (apo)", "Mean (holo)", "Var (apo)", "Var (holo)",
                         "Median (apo)", "Median (holo)", "Last frame (apo)", "Last frame (holo)")



out=paste(
  sim.stats[,1], sim.stats[,2],
  sim.stats[,3], sim.stats[,4],
  sim.stats[,5], sim.stats[,6],
  sim.stats[,7], sim.stats[,8],
  sep=" & ")
out=paste(out, "\\\\")

out=c( paste(paste(colnames(sim.stats), collapse=" & "), "\\\\"), out)

write(out, file="sim.stats.dat")


height <- as.matrix(cbind( sim.stats[,1], sim.stats[,2],
                          sim.stats[,5], sim.stats[,6],
                          sim.stats[,7], sim.stats[,8] ) )


c1="grey60"; c2=1;
col=c(c1,c1,c1,c1,c1,c1,c2,c2,c2,c2,c2,c2)
                         
pdf("sim.stats.pdf", h=4.5,w=6)
mp <- barplot(height, beside=TRUE, col=col, axisnames=FALSE,
              main="Simulation statistics of projections onto PC1")
cm <- colMeans(mp)
at <- c(mean(cm[1:2]), mean(cm[3:4]), mean(cm[5:6]))
mtext(1, at=at, text=c("Mean", "Median", "Last frame"), line=1)
legend(mp[1,1], 250, legend=c("Apo", "Holo"), fill=c(c1, c2))

t1<-t.test(sim.stats[,1], sim.stats[,2], alternative="less")
t2<-t.test(sim.stats[,5], sim.stats[,6], alternative="less")
t3<-t.test(sim.stats[,7], sim.stats[,8], alternative="less")
mtext(1, at=at, text=paste("p=", round(c(t1$p.value,t2$p.value,t3$p.value),3), sep=""), line=2, cex=.7)


dev.off()








hl.col=c("orange", "green", "red")
cex.main=1.2



#pdf("conformerplot_xray_md_ALL_2.pdf", h=6.83, w=6.83)
postscript("conformerplot_xray_md_ALL_2.ps",
           height=10, width=10, paper="special",
           horizontal=FALSE)
par(mfrow=c(4,4), mar=c(2.5,2.5,2,1), mgp=c(1.4,0.5,0) )


conf.plot( d121, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="PC2", xlab="", mtext="A", mtext2="300 ns", main="Closed Unbound", cex.main=1.5 )
conf.plot( d123, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,370), ylim=c(-120,250), ylab="",  xlab="", mtext="B",  mtext2="300 ns", main="Open (r'') Unbound", cex.main=1.5 )
conf.plot( d126, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="",  xlab="", mtext="C",  mtext2="300 ns", main="Closed (t) ATP Bound", cex.main=1.5 )
conf.plot( d124, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,370), ylim=c(-120,250), ylab="", xlab="", mtext="D",  mtext2="300 ns", main="Open ADP (r'') Bound", cex.main=1.5 )


conf.plot( p1, pc.xray, p1a, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="PC2", xlab="" ,mtext="E", mtext2="100 ns", main="Closed (t) Unbound", cex.main=1.5 )
conf.plot( p2,  pc.xray, p1a, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", xlab="", mtext="F",  mtext2="100 ns", main="Closed (t) Unbound", cex.main=1.5 )


conf.plot( p7, pc.xray, p1a, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", xlab="", mtext="K",  mtext2="100 ns", main="Closed (t) ATP Bound", cex.main=1.5 )
conf.plot( p8, pc.xray, p1a, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", xlab="", mtext="L",  mtext2="100 ns", main="Closed (t) ATP Bound", cex.main=1.5 )


conf.plot( p3,  pc.xray, p1a, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="PC2", xlab="PC1", mtext="G",  mtext2="100 ns", main="Closed (t) Unbound", cex.main=1.5 )
conf.plot( p4,  pc.xray, p1a, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", xlab="", mtext="H",  mtext2="100 ns", main="Closed (t) Unbound", cex.main=1.5 )

conf.plot( p9, pc.xray, p1a, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", xlab="PC1", mtext="M",  mtext2="100 ns", main="Closed (t) ATP Bound", cex.main=1.5 )
conf.plot( p10,pc.xray, p1a, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", xlab="", mtext="N",  mtext2="100 ns", main="Closed (t) ATP Bound", cex.main=1.5 )


conf.plot( p5, pc.xray, p1a, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-150,250), ylab="PC2", xlab="PC1", mtext="I", mtext2="100 ns",  main="Closed (t) Unbound", cex.main=1.5 )
conf.plot( p6, pc.xray, p1a, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", xlab="PC1", mtext="J",  mtext2="100 ns", main="Closed (t) Unbound", cex.main=1.5 )


conf.plot( p11,  pc.xray, p1a, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", xlab="PC1", mtext="O",  mtext2="100 ns", main="Closed (t) ATP Bound", cex.main=1.5 )
conf.plot( p12,  pc.xray, p1a, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", xlab="PC1", mtext="P",  mtext2="100 ns", main="Closed(t) ATP Bound", cex.main=1.5 )

dev.off()



pdf("conformerplot_40ns.pdf", h=2)
par(mfrow=c(1,4), mar=c(2.5,2.5,2,1), mgp=c(1.4,0.5,0) )
conf.plot( p5, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", mtext="M", main="(V) Closed Unbound", cex.main=0.8 )
conf.plot( p6, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", mtext="N", main="(VI) Closed Unbound", cex.main=0.8 )
conf.plot( p11, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", mtext="O", main="(V) Closed ATP Bound", cex.main=0.8 )
conf.plot( p12, pc.xray, hl.ind=highlight.ind2, hl.col=hl.col, xlim=c(-150,300), ylim=c(-50,250), ylab="", mtext="P", main="(VI) Closed ATP Bound", cex.main=0.8 )
dev.off()



m1<-c(mean(p1[,1]), mean(p2[,1]), mean(p3[,1]),
      mean(p4[,1]), mean(p5[,1]), mean(p6[,1]) )

m2<-c(mean(p7[,1]), mean(p8[,1]), mean(p9[,1]),
      mean(p10[,1]), mean(p11[,1]), mean(p12[,1]) )

t.test(m1, m2, var.equal=T, alternative="less")


var1<-c(var(p1[,1]), var(p2[,1]), var(p3[,1]),
      var(p4[,1]), var(p5[,1]), var(p6[,1]) )

var2<-c(var(p7[,1]), var(p8[,1]), var(p9[,1]),
      var(p10[,1]), var(p11[,1]), var(p12[,1]) )







conf.plot <- function( projection, pc.xray, projection.secondary=NULL,
                      pc.inds = c(1,2),
                      xlim=NULL, ylim=NULL, main = "", xlab="PC1", ylab="PC2",
                      hl.col=NULL, hl.ind=NULL, ablines=NULL,
                      mtext=NULL, mtext2=NULL, ... ) {

  
  if( is.null(projection.secondary) ) {
    projection.secondary = projection
    col.sec = densCols(projection[,pc.inds[1]], projection[,pc.inds[2]])
    projection = NULL
  }
  else {
    col.sec="grey90"
  }
    
  plot(projection.secondary[,pc.inds[1]], projection.secondary[,pc.inds[2]], type="p",
       main = main, xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim,
       col=col.sec,
       pch=20, cex.axis=1.24, ...)

  if( !is.null(projection) ) {
    points( projection[,pc.inds[1]], projection[,pc.inds[2]],
           col=densCols(projection[,pc.inds[1]], projection[,pc.inds[2]]),
           pch=20)
  }
    
  points( pc.xray$z[,pc.inds[1]], pc.xray$z[,pc.inds[2]], col = "black", pch = 16, cex = 1 )
  points( pc.xray$z[,pc.inds[1]], pc.xray$z[,pc.inds[2]], col="gray", pch = 16, cex = 0.4 )
  points( pc.xray$z[hl.ind, pc.inds[1]], pc.xray$z[hl.ind, pc.inds[2]], col = hl.col, pch = 16 )
  
  mtext(mtext, side=3, at=xlim[1]-50, cex=1., line=0.6)
  #mtext(mtext2, side=3, at=xlim[1]*1.2, cex=0.6, line=-1)
  text(xlim[1], ylim[2]*0.95, mtext2, pos=4, cex=1.1)

  #abline(v=mean(projection.secondary[,pc.inds[1]]), col="grey80", lty=2)

}