library(bio3d)
library(car)

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




pdf("xray_confplot_EQ_And_Tstates.pdf", w=12, h=6)



load("xrayPCA_Tstates.RData")
id <- substr(basename(pdbs$id), 1, 6)


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_M", id),
                    
                    grep("1SX3_A", id),grep("1KP8_B", id),
                    grep("1KP8_F", id),grep("1SX3_G", id),
                    grep("1SX3_C", id),grep("2C7C_H", 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
  )

t1 <- which(fit$cluster==1)
t2 <- which(fit$cluster==2)
t3 <- which(fit$cluster==3)
cluster.col=fit$cluster
cluster.col[t1]="orange"
cluster.col[t2]="red"
cluster.col[t3]="green"

cex.axis=1.1
cex.lab=1.1
cex=0.9

par(mgp=c(2,1,0), mfcol=c(1,2))
plot(pc.xray$z[,1], pc.xray$z[,2], col=cluster.col, xlab = "PC1", ylab = "PC2", 
        ylim=c(-35,23), xlim=c(-40,50), cex=cex, pch=20,
     cex.axis=cex.axis, cex.lab=cex.lab)
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.5, 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")
mtext("A", at=-40, cex=1.75)

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)





load(file="../13_xrayPCA_EQ/xrayPCA_EQ.RData")
id <- substr(basename(pdbs$id), 1, 6)


highlight.ind <-  c(
                    grep("1GRU_M", id), grep("1SX4_N", id),
                    grep("1SVT_H", id),
                         
                    grep("1KP8_B", id),grep("1SS8_B", id),
                    grep("1XCK_B", id),grep("3FBH_A", id),
                    grep("2NWC_I", id),

                    grep("1SX4_A", id), grep("1SVT_D", id), 
                    grep("1AON_A", id),

                    grep("1SX3_B", id), grep("1SX3_N", id),
                    grep("3E76_H", id), 

                    grep("2C7E_A", id),

                    grep("2C7C_K", id),
                    grep("2C7C_J", id), grep("2C7C_H", id),
                    grep("2C7C_L", id), 
                    grep("2C7C_B", id), 
                    grep("2C7C_A", id), 
                    
                    grep("2C7D_I", id),
                     grep("2C7D_A", id),
                    grep("2C7D_L", id), grep("2C7D_M", id) )


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

fit<-kmeans( cbind(pc.xray$z[,1], pc.xray$z[,2]), 4, 2)

t1 <- which(fit$cluster==1)
t2 <- which(fit$cluster==2)
t3 <- which(fit$cluster==3)
t4 <- which(fit$cluster==4)
cluster.col=fit$cluster
cluster.col[t1]="red"
cluster.col[t2]="orange"
cluster.col[t3]="green"
cluster.col[t4]="yellow"


plot(pc.xray$z[,1], pc.xray$z[,2], col=cluster.col, xlab = "PC1", ylab = "PC2", 
	ylim=c(-12,12), xlim=c(-25,27), cex=cex, cex.axis=cex.axis, cex.lab=cex.lab, pch=20)
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.5, 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")
mtext("B", at=-25, cex=1.75)

radius=c(4,4,3,3,2)

## T conf
for ( i in 1:4 ) {
  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(1, 9, "GroEL (apo)", cex=.9)
text(-6.5,-9.5, "GroEL-ATP (cis)", cex=.9)
text(-17,6, "GroEL-ADP (cis)", cex=.9)
text(18,1.5, "GroEL-ADP (trans)", cex=.9)
dev.off()