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

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


##### 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"
anchor.res = c(51, 243:248)  ## From the core calculation

setwd(output_dir)

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)
pdb.names = gsub('.{4}$', '', unlist(pdb.namelist))  

## align to the anchor residues calculated from Bio3D Core Calculation
pdb = read.pdb(pdb_ref)
ca.inds = atom.select(pdb, resno = anchor.res)
## align against the PDB and anchor residues
xyz <- fit.xyz(fixed = pdb$xyz, mobile=pdb.df,
               fixed.inds=ca.inds$xyz,
               mobile.inds=ca.inds$xyz)
dim(xyz)

#### 1. Plot RMSF #######

resno = pdb$atom$resno
rf = rmsf(xyz, grpby=resno)
pdf(file = paste(clust.type, "_rmsf", ".pdf", sep=""))
plot.bio3d(rf, ylab="RMSF (Å)", xlab="Residue No.", typ="l", ylim = c(0,60))
points(rf)
dev.off()
rf.df = data.frame(res = 1:length(rf), rmsf = rf)
write.csv(rf.df, file = paste(clust.type, "_rmsf", ".csv", sep=""), row.names = FALSE)


##### 2. PCA clustering  ###########

pc.xray = pca.xyz(xyz, use.svd = TRUE)

## structure wise PC contributions
pdf(file = paste(clust.type, "_pca_plot.pdf", sep=""))
plot(pc.xray)
dev.off()


hc.pca <- hclust(dist(pc.xray$z[,1:5]))
#hc.pca <- hclust(dist(pc.xray$z[,1:2]))

# Dendrogram plot
# names(cols) <- pdb.names
pdf(file = paste(clust.type, "_pca_clust.pdf", sep=""))
hclustplot(hc.pca, k=4, labels = pdb.names, ylab="Distance in PC Space", main="PC1-5",
           cex=0.37, fillbox=FALSE)
dev.off()

## identify the clusters
group.pclabels = cutree(hc.pca, k = 4)
pc.clust1 = pdb.names[group.pclabels == 1]
pc.clust2 = pdb.names[group.pclabels == 2]
pc.clust3 = pdb.names[group.pclabels == 3]
pc.clust4 = pdb.names[group.pclabels == 4]
length(pc.clust1); length(pc.clust2); length(pc.clust3); length(pc.clust4)
save(pc.clust1, pc.clust2, pc.clust3, pc.clust4, 
     file = "all_pca_clust_groups.Rda")

## 3. Structure wise PC components and 3D plot of PCA components
pc.df = as.data.frame(pc.xray$z)
colnames(pc.df) = paste(rep("PC", ncol(pc.df)), 1:ncol(pc.df), sep = "")
pc.df$pdbnames = pdb.names
pc.df$class = as.character(group.pclabels)
dim(pc.df); class(pc.df); class(pc.df$class)
## pc.df contains residue wise PC contribution, rows are PDB structures, and columns are PC components
write.table(pc.df, file = paste(clust.type, "_PCcomponents.txt"), sep="\t",
            quote = FALSE, row.names = FALSE)
# plot(pc.df$PC1, pc.df$PC2)

library(plotly)
pc123_3D = plot_ly(pc.df, x = ~PC1, y=~PC2, z=~PC3, color = ~class, mode = "markers")
htmlwidgets::saveWidget(as_widget(pc123_3D), "pc123_3D_plot.html")
pc124_3D =  plot_ly(pc.df, x = ~PC1, y=~PC2, z=~PC4, color = ~class)
htmlwidgets::saveWidget(as_widget(pc124_3D), "pc124_3D_plot.html")


#### 4. Per residue contribution to PCA   ###########3

## Find normalized contribution of x,y,z coordiantes from pc.xray$U
## normalized eigenvectors = sqrt(x^2 + y^2 + z^2)
xyz.eigenvec = pc.xray$U
seq.vec = 1:nrow(xyz.eigenvec)
seq.list = split(seq.vec, ceiling(seq_along(seq.vec)/3))
normalize = function(y){
  #print(xyz.eigenvec[y, 1:5])
  return(apply(xyz.eigenvec[y, ], 2, FUN = function(x) sqrt(sum(x^2))))
}
norm.eigen = lapply(seq.list, normalize)
norm.eigen.mat = do.call(rbind, norm.eigen) 
dim(norm.eigen.mat)

### Average atom-wise PC contributions to per residue contributions
res.wt = list()
for (resnum in unique(pdb$atom$resno)){
  ca.inds = atom.select(pdb, resno = resnum)
  res.wt[[resnum]] = apply(norm.eigen.mat[ca.inds$atom, ], 2, mean)
}
res.wt.mat = do.call(rbind, res.wt) 
dim(res.wt.mat)

res.wt.mat = as.data.frame(res.wt.mat)
colnames(res.wt.mat) = paste(rep("PC", ncol(res.wt.mat)), 1:ncol(res.wt.mat), sep = "")
res.wt.mat$residue = 1:nrow(res.wt.mat)
## move residue column to first
# res.wt.mat <- res.wt.mat[,c(ncol(res.wt.mat),1:ncol(res.wt.mat)-1)]
write.table(res.wt.mat, file = paste(clust.type, "RESIDUEwise_PCA_Contribution.txt", sep = "_"), sep="\t",
            quote = FALSE, row.names = FALSE)

pdf(file = paste(clust.type, "RESIDUEwise_PC_Contribution.pdf", sep="_"))
par(mfrow = c(3, 1),cex = 0.75, mar = c(3, 4, 1, 1))
plot(1:nrow(res.wt.mat), res.wt.mat[,1], ylim = c(0.000, 0.050), type = "l", col = "blue", ylab = "PC1")
plot(1:nrow(res.wt.mat), res.wt.mat[,2], ylim = c(0.000, 0.050), type = "l", col = "blue", ylab = "PC2")
plot(1:nrow(res.wt.mat), res.wt.mat[,3], ylim = c(0.000, 0.050), type = "l", col = "blue", ylab = "PC3")

par(mfrow = c(3, 1),cex = 0.75, mar = c(3, 4, 1, 1))
plot(1:nrow(res.wt.mat), res.wt.mat[,4], ylim = c(0.000, 0.050), type = "l", col = "blue", ylab = "PC4")
plot(1:nrow(res.wt.mat), res.wt.mat[,5], ylim = c(0.000, 0.050), type = "l", col = "blue", ylab = "PC5")
plot(1:nrow(res.wt.mat), res.wt.mat[,6], ylim = c(0.000, 0.050), type = "l", col = "blue", ylab = "PC6")
dev.off()



## 5. Write out interpolated trajectory for various PC components only considering backbone atoms  #####
ca.inds.bb = atom.select(pdb, elety=c("C3'", "C4'", "C5'", "P"))
# xyz <- fit.xyz(fixed = pdb$xyz, mobile=pdb.df,
#                fixed.inds=ca.inds.bb$xyz,
#                mobile.inds=ca.inds.bb$xyz)
# dim(xyz)
pc.xray = pca.xyz(xyz[, ca.inds.bb$xyz], use.svd = TRUE)
mktrj.pca(pc.xray, pc=1, pdb = trim(pdb, ca.inds.bb), step = 0.3, file = "pc1_bb_traj.pdb")
mktrj.pca(pc.xray, pc=2, pdb = trim(pdb, ca.inds.bb), step = 0.3, file = "pc2_bb_traj.pdb")
mktrj.pca(pc.xray, pc=3, pdb = trim(pdb, ca.inds.bb), step = 0.3, file = "pc3_bb_traj.pdb")
mktrj.pca(pc.xray, pc=4, pdb = trim(pdb, ca.inds.bb), step = 0.3, file = "pc4_bb_traj.pdb")
mktrj.pca(pc.xray, pc=5, pdb = trim(pdb, ca.inds.bb), step = 0.3, file = "pc5_bb_traj.pdb")







