library(bio3d)

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

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

hits <- plot.blast(blast, cutoff=500 ) 

##raw.files <- get.pdb(unq.ids, path = "raw_pdbs_new")
unq.ids <- unique(substr(hits$pdb.id, 1, 4))
files <- paste("raw_pdbs_new/", unq.ids, ".pdb", sep = "")

library(multicore)
threads <- 6
ptm <- proc.time()
files.multicore <- rep(1:threads, each=length(files)/threads)


files.split <- c()
path <- "raw_pdbs_new/split_chain_dimers"
jobs <- list()

for ( i in 1:threads ) {
  job.inds <- which(files.multicore == i)
  q <- parallel(start.split(files[job.inds], path))
  jobs[[i]]=q
}

res=collect(jobs, wait=TRUE)
names <- c()
for (job in res) {
  names=c(names, job$names)
}

files.inds = grep("1AON_[HIJKLMN]|1SX4_[HIJKLMN]|1PF9_[HIJKLMN]|1SVT_[HIJKLMN]|1XCK|1KP8|1SX3|2NWC_[HIJKLMN]|3FBH|2C7E_[ABCDEFG]", files.split)

files.inds = grep("1AON_[HIJKLMN]|1SX4_[HIJKLMN]|1PF9_[HIJKLMN]|1SVT_[HIJKLMN]|1XCK|1KP8|1SX3|2NWC_[HIJKLMN]|3FBH", files.split)

## perform alignment
aln <- pdbaln(files.split[files.inds])

## Core superposition ##
pdbs <- read.fasta.pdb(aln, "", "", het2atom = T)
core <- core.find(pdbs)

print(core, 1)



##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/", full.pdbs = TRUE, het2atom = TRUE)


load(file="monomer.core.RData")
core.md <- c(monomer.core$c1A.xyz, monomer.core$c1A.xyz+1572)
core.exp <- c(core$c1A.xyz, core$c1A.xyz+1575)


#xyz <- fit.xyz(fixed = pdbs$xyz[1,], mobile = pdbs, fixed.inds = test.core , mobile.inds = test.core, pdb.path = "", pdbext = "", outpath = "core_fitlsq/", full.pdbs = TRUE, het2atom = TRUE)
xyz <- fit.xyz(fixed = pdbs$xyz[1,], mobile = pdbs, fixed.inds = core.exp , mobile.inds = core.exp, full.pdbs = FALSE, het2atom = TRUE)


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


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

pdf("pc.pdf")
plot(pc.xray)
dev.off()



save(pdb.reference, blast, files, hits, unq.ids, files.split, aln, pdbs, xyz,
     files.inds, gaps.pos, gaps.res, pc.xray, monomer.core
     core.exp, core.md, file="pca_dimers.RData")

save(pdb.reference, blast, files, hits, unq.ids, files.split, aln, pdbs, xyz,
     files.inds, gaps.pos, gaps.res, pc.xray, monomer.core,
     core.exp, core.md, file="pca_dimers_t-only.RData")




##mag=2, step=0.25,
a <- mktrj.pca(pc.xray, pc=1, file="pc1_tr.pdb", 
               resno = pdbs$resno[1, gaps.res$f.inds],
               resid = aa123(pdbs$ali[1, gaps.res$f.inds]) )
b <- mktrj.pca(pc.xray, pc=2, file="pc2_tr.pdb",mag=2, step=0.25,
               resno = pdbs$resno[1, gaps.res$f.inds],
               resid = aa123(pdbs$ali[1, gaps.res$f.inds]) )
a <- mktrj.pca(pc.xray, pc=1, file="pc1_t.pdb", 
               resno = pdbs$resno[1, gaps.res$f.inds],
               resid = aa123(pdbs$ali[1, gaps.res$f.inds]) )
b <- mktrj.pca(pc.xray, pc=2, file="pc2_t.pdb",mag=2, step=0.25,
               resno = pdbs$resno[1, gaps.res$f.inds],
               resid = aa123(pdbs$ali[1, gaps.res$f.inds]) )




b <- mktrj.pca(pc.xray, pc=3, file="pc3.pdb",
               resno = pdbs$resno[1, gaps.res$f.inds],
               resid = aa123(pdbs$ali[1, gaps.res$f.inds]) )

b <- mktrj.pca(pc.xray, pc=4, file="pc4.pdb",
               resno = pdbs$resno[1, gaps.res$f.inds],
               resid = aa123(pdbs$ali[1, gaps.res$f.inds]) )



id <- substr(basename(pdbs$id), 1, 6)


pdf("xray_confplot_t-only.pdf", w=6, h=6)
par(mgp=c(2,1,0))
plot(pc.xray$z[,1], pc.xray$z[,2], xlab = "PC1", ylab = "PC2", 
        cex=0.8, pch=20, cex.axis=.8, cex.lab=.8)
text(pc.xray$z[, 1], pc.xray$z[, 2], id, cex=0.6,  col="grey50")
##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")
abline(h=0, lty=2, col="grey60")
abline(v=0, lty=2, col="grey60")

dev.off()













if(FALSE){
files.split <- c()
path <- "raw_pdbs_new/split_chain_dimers"
for ( i in 1:30 ) {
  pdb <- read.pdb(files[i], maxlines=99999)
  chains <- unique(pdb$atom[, "chain"])

  for ( j in 1:6 ) {
    chain.inds <- c(j, j+1)
    new.chains <- chains[chain.inds]
    sel <- atom.select(pdb, chain=new.chains)
    new.pdb <- trim.pdb(pdb, sel)
    new.name <- paste(substr(basename(files[i]), 
                             1, 4), "_", paste(new.chains, collapse=""), ".pdb", sep = "")
    new.name <- file.path(path, new.name)
    write.pdb(new.pdb, file = new.name)
    files.split <- c(files.split, new.name)
    }

  if ( "H" %in% chains ) {
    for ( j in 8:13 ) {
      chain.inds <- c(j, j+1)
      new.chains <- chains[chain.inds]
      sel <- atom.select(pdb, chain=new.chains)
      new.pdb <- trim.pdb(pdb, sel)
      new.name <- paste(substr(basename(files[i]), 
                               1, 4), "_", paste(new.chains, collapse=""), ".pdb", sep = "")
      new.name <- file.path(path, new.name)
      write.pdb(new.pdb, file = new.name)
      files.split <- c(files.split, new.name)
    }
  }
}
}






if(FALSE){
i <- grep("2C7C", files)
pdb <- read.pdb(files[i], maxlines=99999)
chains <- unique(pdb$atom[, "chain"])

for ( j in 1:6 ) {
  chain.inds <- c(j, j+1)
  new.chains <- chains[chain.inds]

  sel <- atom.select(pdb, chain=new.chains)
  ref.pdb <- trim.pdb(pdb, sel)
  
  sel.a <- atom.select(pdb, chain=new.chains[1])
  sel.b <- atom.select(pdb, chain=new.chains[2])
  pdb.a <- trim.pdb(pdb, sel.a)
  pdb.b <- trim.pdb(pdb, sel.b)

  #db.a$atom[,"chain"] = rep(new.chains[2], length(pdb.a$atom[,"chain"]))
  #db.b$atom[,"chain"] = rep(new.chains[1], length(pdb.b$atom[,"chain"]))
  

  new.pdb <- ref.pdb
  #new.pdb$atom <- rbind(pdb.b$atom, pdb.a$atom)
  #new.pdb$xyz <- rbind(pdb.b$xyz, pdb.a$xyz)
  #new.pdb$atom[,"eleno"] <- seq(1, nrow(new.pdb$atom))
  #new.pdb$atom[,"chain"] <- c(pdb.a$atom[,"chain"], pdb.b$atom[,"chain"])
  
  new.name <- paste(substr(basename(files[i]), 
                           1, 4), "_", paste(new.chains, collapse=""), ".pdb", sep = "")
  
  new.name <- file.path(path, new.name)
  write.pdb(new.pdb, xyz=as.numeric(c(pdb.a$xyz, pdb.b$xyz)),
            end=FALSE, file=new.name)
            ##chain=c(pdb.a$atom[,"chain"], pdb.b$atom[,"chain"]),
            ##end=FALSE, ##eleno=new.pdb$atom[,"eleno"])
  
  print(new.name)
  ##files.split <- c(files.split, new.name)
}
}