## Analyze the top PDB sturctures throught RMSD and PCA clustering

# install.packages("stringr")
# install.packages("plyr")
library(bio3d)
library(stringr)
library(tidyr)
library(plyr)
library(dplyr)


##### Begin Analysis #######

### INPUTS to provide: output prefix (clust.type), output_dir, pdb_path(path), pdb_to_align (pdb_ref)

##### Read through the PDBs and store them in the XYZ matrix (pdb.df) for further analysis ######
clust.type = "all"
output_dir = "/Users/bhandariyr/Documents/projects/afm/rnaseP/PCA/FinalAnalysis/July2021_set2/overall"
path = "/Users/bhandariyr/Documents/projects/afm/rnaseP/PCA/FinalAnalysis/July2021_set2/data/FMT_PDB" 
pdb_ref = "/Users/bhandariyr/Documents/projects/afm/rnaseP/PCA/FinalAnalysis/July2021_set2/data/FMT_PDB/9758_P78_fmt.pdb"

setwd(output_dir)  ## set working directory
getwd()

## Sanity Check: compare exact atom numbering between crystal PDB and cafemol fitted allatom PDBs
pdb1 = read.pdb(pdb_ref)
pdb2 = read.pdb("/Users/bhandariyr/Documents/projects/afm/rnaseP/PCA/FinalAnalysis/July2021_set2/data/FMT_PDB//RNaseP_tr_fmt.pdb")

## check if atom numbering are identical, all lines below should return TRUE
identical(pdb1$atom$elety, pdb2$atom$elety)
identical(pdb1$atom$resid, pdb2$atom$resid)
identical(pdb1$atom$resno, pdb2$atom$resno)

## Read PDB coordinates into matrix  #######

files = list.files(path = path, pattern = ".pdb" )
pdb.namelist = list()
pdb.dflist = list()
count = 1 ## counter for list
for (file in files){
  file.path = paste(path, file, sep = "/")
  pdb = read.pdb(file.path)
  pdb.namelist[[count]] = file
  pdb.dflist[[count]] = pdb$xyz[1, ]
  count = count + 1
  #print((pdb$xyz[1, ]))
}

pdb.df = do.call(rbind, pdb.dflist) 
dim(pdb.df)
# View(pdb.df)
pdb.names = gsub('.{4}$', '', unlist(pdb.namelist))

## alignment and rmsd clustering
#
pdb = read.pdb(pdb_ref)
# identical(pdb$atom$elety, pdb2$atom$elety)
## remove the missing residues in the crystal structure
missing.res = c(82:100, 115:121, 203:226, 266:268)
xtal.res = setdiff(c(1:268), missing.res)
ca.inds = atom.select(pdb, resno = xtal.res)
# xtal.res = c(1:10, 10:20) ## custom residues
pdb$xyz
## align against the given reference PDB
xyz <- fit.xyz(fixed = pdb$xyz, mobile=pdb.df,
               fixed.inds=ca.inds$xyz,
               mobile.inds=ca.inds$xyz)



#### Core Residue Calculation  ####

## select only "P" atoms in the XYZ, add coordinates 1:3 for first residue
ca.inds = atom.select(pdb, elety = "P") 
pdb.plist = list()
count2 = 1 ## counter for list
for (modelNum in 1:nrow(xyz)){
  pdb.plist[[count2]] = xyz[modelNum, ][c(1:3, ca.inds$xyz)]
  count2 = count2 + 1
}
pdb.df.P = do.call(rbind, pdb.plist) 
dim(pdb.df.P)

core = core.find(pdb.df.P, stop.at = 2)

pdf(file = "CoreVolume.pdf")
col=rep("black", length(core$volume))
col[core$volume<50]="pink"; col[core$volume<1]="red"
# plot(core$length, core$volume, col=col)
plot(core, col = col)
#plot(core$length, core$volume, xlim=c(268, 0))
dev.off()

save(core, file = paste(clust.type, "CoreVolume.Rda", sep = "_"))
core.inds <- print(core, vol=1)  ## residues with volume less than 1 A3
core.inds = print(core, vol = 50)
# # 11 positions (cumulative volume <= 50 Angstrom^3) 
# start end length
# 1    49  51      3
# 2   242 249      8


core.df = data.frame(core_size = core$length, volume = core$volume)
write.csv(core.df, file = paste(clust.type, "core_paramters", ".csv", sep="_"), row.names = FALSE)







