library(bio3d)
##load("xrayPCA.RData")
source("confplot_funs.R")
library(car)

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 )

unq.ids <- unique(substr(hits$pdb.id, 1, 4))
files <- paste("raw_pdbs_new/split_chain/", hits$pdb.id, ".pdb", sep = "")
exclude.inds = grep("1AON_[HIJKLMN]|1SX4_[HIJKLMN]|1PF9_[HIJKLMN]|1SVT_[HIJKLMN]|1XCK|1KP8|1SX3|1SS8|2NWC_[HIJKLMN]|3FBH|2C7E_[ABCDEFG]", hits$pdb.id)
files.excl <- files[exclude.inds]

##2C7C_[HIJKLMN]|1GRU_[HIJKLMN]

aln <- pdbaln(files.excl)
pdbs <- read.fasta.pdb(aln, "", "", het2atom = T)

core <- core.find(pdbs)

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

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

pc.xray <- pca.xyz(xyz[,gaps.pos$f.inds])

save(pdb, hits, unq.ids, files, files.excl, aln, pdbs, core, gaps.pos, gaps.res, xyz, pc.xray, file="xrayPCA.RData")


load("xrayPCA.RData")
vi <- mktrj.pca(pc.xray, pc=1, file="pc1_eq.pdb")
vi <- mktrj.pca(pc.xray, pc=2, file="pc2_eq.pdb", mag=4, step=0.50)
vi <- mktrj.pca(pc.xray, pc=3, file="pc3_eq.pdb", mag=4, step=0.50)




library(gclus)
library(maptools)

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)
fit<-kmeans( cbind(pc.xray$z[,1], pc.xray$z[,2]), 3, 20)


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


highlight.ind <-  c(grep("1SX4_N", id),grep("1PF9_N", id),
                    grep("1AON_K", id),grep("1SVT_I", id),
                    grep("1SVT_K", id),grep("1AON_J", id),
                    
                    grep("1SX3_A", id),grep("1KP8_B", id),
                    grep("1KP8_F", id),grep("1SX3_I", id),
                    grep("1SX3_C", id),grep("2C7C_H", id),

                    grep("2C7E_A", id),grep("2C7E_G", id),

                    grep("1SS8_C", id),grep("3FBH_N", id),
                    grep("1XCK_B", id),grep("1XCK_D", id),
                    grep("1SS8_A", id),grep("1XCK_N", id) )


                    
pos=c(
  4,4,4,4,4,1,
  4,4,4,1,4,4,
  4,4,
  4,4,4,4,4,4
  )


pdf("xray_confplot.pdf", w=6, h=12)
par(mfcol=c(2,1), mgp=c(2,1,0))
plot(pc.xray$z[,1], pc.xray$z[,2], col=fit$cluster, xlab = "PC1", ylab = "PC2", 
        ylim=c(-30,40), xlim=c(-50,200), cex=0.8, pch=20, cex.axis=.8, cex.lab=.8)
##text(pc.xray$z[highlight.ind, 1], pc.xray$z[highlight.ind, 2], id[highlight.ind], cex=0.6, pos=pos, col="grey50")
text(pc.xray$z[, 1], pc.xray$z[, 2], id, cex=0.6, pos=pos, col="grey50")
abline(h=0, lty=2, col="grey60")
abline(v=0, lty=2, col="grey60")

radius=c(2,2,2,2,2)

## T conf
for ( i in 1:3 ) {
  r2.inds <- which(fit$cluster==i)
  #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), radius[i], lty="dashed", col="grey80", center.pch=0, lwd=0.7)
  #text(mean(v1)-5, mean(v2)+5, "Cis/Trans t", cex=0.9, col="black")
}

text(-8, -10, "GroEL (apo)", cex=.9)
text(18,12, "GroEL-ATP (cis)", cex=.9)
text(-25,7, "GroEL-ADP (trans)", cex=.9)

plot(pc.xray$z[,3], pc.xray$z[,2], col=fit$cluster, xlab = "PC3", ylab = "PC2", 
        ylim=c(-30,40), xlim=c(-20,40), cex=0.8, pch=20, cex.axis=.8, cex.lab=.8)
##text(pc.xray$z[highlight.ind, 1], pc.xray$z[highlight.ind, 2], id[highlight.ind], cex=0.6, pos=pos, col="grey50")
text(pc.xray$z[, 3], pc.xray$z[, 2], id, cex=0.6, pos=pos, col="grey50")
abline(h=0, lty=2, col="grey60")
abline(v=0, lty=2, col="grey60")

radius=c(2,2,2,2,2)

dev.off()