## Find the top Correlation-Energy (ce) Score structures for the top kappa cohort
## Usage: Rscript input_folder top_percent pdbfile rmsd_file 

library(dplyr)
library(stringr)
library(bio3d)
#install.packages("gridExtra")
library(gridExtra)

args = commandArgs(trailingOnly=TRUE)
#print(args[1])
#print(length(args))
if (length(args) != 4) {
  stop("Please provide the input_folder top_percent pdbfile rmsd_file", call.=FALSE)
}

#cohort = "/Users/bhandariyr/Documents/projects/afm/cobalamin/miscen/scripts/stats/contact10_5/output1_0"
cohort = args[1]
top_percent = as.numeric(args[2])
pdbfile = args[3]
rmsd_path = args[4]

ext1 = paste(cohort,"*.ts", sep = "/")
file.ts = Sys.glob(ext1)
ext2 = paste("top_",top_percent, "%_pdbs", sep = "")
#rmsd_path = paste(cohort,ext2,"rmsd.txt", sep = "/")
#print(rmsd_path)
df.rmsd = read.table(rmsd_path, header = TRUE)
#head(df.rmsd)
## Read and clean timestep file
data.ts.lines = readLines(file.ts)[grep("#all", readLines(file.ts))]
data.ts =  data.ts.lines %>%
  strsplit("\\s+") %>%
  do.call(rbind, .) %>%
  data.frame()
dim(data.ts)
colnames(data.ts) = unlist(strsplit(readLines(file.ts)[grep("#unit", readLines(file.ts))], "\\s+"))
#head(data.ts)
colnames(data.ts)[12] = "stack"
# View(data.ts)
data.ts$frame = as.numeric(rownames(data.ts)) - 1

data.ts.sort <- data.ts %>% select("frame", everything())
## subset to remove frame zero
data.ts.sort = data.ts.sort[2:nrow(data.ts.sort),]
## add the rmsd column
data.ts.sort$calc_rmsd = df.rmsd$rmsd

## Calculate the correlation energy (ce) score, and add to the data frame
afmcc = as.numeric(data.ts.sort$afmcc) #Afmcc score
#toten = as.numeric(data.ts.sort$etot)  #total energy
toten = as.numeric(data.ts.sort$local)+as.numeric(data.ts.sort$go)+as.numeric(data.ts.sort$repul)+
           as.numeric(data.ts.sort$stack)+as.numeric(data.ts.sort$hbond)+as.numeric(data.ts.sort$elect)

cc_norm = (afmcc - min(afmcc))/(max(afmcc) - min(afmcc))
toten_norm = (toten - max(toten))/(max(toten) - min(toten))
data.ts.sort$ce_score = cc_norm*toten_norm
#View(data.ts.sort)

toten_local = as.numeric(data.ts.sort$local)
toten_go = as.numeric(data.ts.sort$go)
toten_repul = as.numeric(data.ts.sort$repul)
toten_stack = as.numeric(data.ts.sort$stack)
toten_hbond = as.numeric(data.ts.sort$hbond)
toten_elect = as.numeric(data.ts.sort$elect)
toten_afmfit = as.numeric(data.ts.sort$afmfit)
toten_stage = as.numeric(data.ts.sort$stage)

toten_norm_local = (toten_local - max(toten_local))/(max(toten_local) - min(toten_local))
toten_norm_go = (toten_go - max(toten_go))/(max(toten_go) - min(toten_go))
toten_norm_repul = (toten_repul - max(toten_repul))/(max(toten_repul) - min(toten_repul))
toten_norm_stack = (toten_stack - max(toten_stack))/(max(toten_stack) - min(toten_stack))
toten_norm_hbond = (toten_hbond - max(toten_hbond))/(max(toten_hbond) - min(toten_hbond))
toten_norm_elect = (toten_elect - max(toten_elect))/(max(toten_elect) - min(toten_elect))
toten_norm_afmfit = (toten_afmfit - max(toten_afmfit))/(max(toten_afmfit) - min(toten_afmfit))
toten_norm_stage = (toten_stage - max(toten_stage))/(max(toten_stage) - min(toten_stage))

data.ts.sort$ce_local = cc_norm*toten_norm_local
data.ts.sort$ce_go = cc_norm*toten_norm_go
data.ts.sort$ce_repul = cc_norm*toten_norm_repul
data.ts.sort$ce_stack = cc_norm*toten_norm_stack
data.ts.sort$ce_hbond = cc_norm*toten_norm_hbond
data.ts.sort$ce_elect = cc_norm*toten_norm_elect
data.ts.sort$ce_afmfit = cc_norm*toten_norm_afmfit
data.ts.sort$ce_stage = cc_norm*toten_norm_stage

## Remove the first frame as it is almost an outlier
data.ts.sort = data.ts.sort[-1,]

## Get the top X percent structures based on ce score
top.frames = round(top_percent*nrow(data.ts.sort)/100)
data.ts.rank = data.ts.sort[order(data.ts.sort$ce_score, decreasing = FALSE), ]
data.ts.top = data.ts.rank[1:top.frames, ]
#View(data.ts.top)

## Do the plotting
dir_name = paste("top_", top_percent, "percent_ce_newEnergy_pdbs", sep = "")
dir.create(file.path(dir_name), recursive = FALSE, showWarnings = TRUE)

plot.name = paste("top",top_percent,"percent_ce_score_newEn_plots.pdf",sep = "")
plot.path = paste(dir_name, plot.name, sep = "/")
pdf(plot.path)
library(ggplot2)
p1 = ggplot(data.ts.sort, aes(x=ce_score, y=calc_rmsd)) +
  geom_point() +
  theme_classic(base_size = 10)+ ylab("calc_rmsd") + 
  ggtitle("all frames") +
  xlab("ce_score")

p2 = ggplot(data.ts.top, aes(x=ce_score, y=calc_rmsd)) +
  geom_point() +
  theme_classic(base_size = 10)+ ylab("calc_rmsd") + 
  ggtitle(paste("top_", top_percent, "_percent_frames", sep = "")) +
  xlab("ce_score")

data.ts.sort$afmcc = as.numeric(data.ts.sort$afmcc)
data.ts.top$afmcc = as.numeric(data.ts.top$afmcc)

p3 = ggplot(data.ts.sort, aes(x=afmcc, y=calc_rmsd)) +
  geom_point() +
  theme_classic(base_size = 10)+ ylab("calc_rmsd") +
  ggtitle("all frames") +
  xlab("afmcc")

data.ts.sort.afmcc = data.ts.sort[order(data.ts.sort$afmcc, decreasing = TRUE), ]
data.ts.top.afmcc = data.ts.sort.afmcc[1:60,]
p4 = ggplot(data.ts.top.afmcc, aes(x=afmcc, y=calc_rmsd)) +
  geom_point() +
  theme_classic(base_size = 10)+ ylab("calc_rmsd") +
  ggtitle(paste("top_", top_percent, "_percent_frames", sep = "")) +
  xlab("afmcc")

p5 = ggplot(data.ts.sort, aes(x=ce_local, y=calc_rmsd)) +
  geom_point() +
  theme_classic(base_size = 10)+ ylab("calc_rmsd") +
  ggtitle("all frames") +
  xlab("ce_local")

p6 = ggplot(data.ts.sort, aes(x=ce_go, y=calc_rmsd)) +
  geom_point() +
  theme_classic(base_size = 10)+ ylab("calc_rmsd") +
  ggtitle("all frames") +
  xlab("ce_go")

p7 = ggplot(data.ts.sort, aes(x=ce_repul, y=calc_rmsd)) +
  geom_point() +
  theme_classic(base_size = 10)+ ylab("calc_rmsd") +
  ggtitle("all frames") +
  xlab("ce_repul")

p8 = ggplot(data.ts.sort, aes(x=ce_stack, y=calc_rmsd)) +
  geom_point() +
  theme_classic(base_size = 10)+ ylab("calc_rmsd") +
  ggtitle("all frames") +
  xlab("ce_stack")

p9 = ggplot(data.ts.sort, aes(x=ce_hbond, y=calc_rmsd)) +
  geom_point() +
  theme_classic(base_size = 10)+ ylab("calc_rmsd") +
  ggtitle("all frames") +
  xlab("ce_hbond")

p10 = ggplot(data.ts.sort, aes(x=ce_elect, y=calc_rmsd)) +
  geom_point() +
  theme_classic(base_size = 10)+ ylab("calc_rmsd") +
  ggtitle("all frames") +
  xlab("ce_elect")

p11 = ggplot(data.ts.sort, aes(x=ce_afmfit, y=calc_rmsd)) +
  geom_point() +
  theme_classic(base_size = 10)+ ylab("calc_rmsd") +
  ggtitle("all frames") +
  xlab("ce_afmfit")

p12 = ggplot(data.ts.sort, aes(x=ce_stage, y=calc_rmsd)) +
  geom_point() +
  theme_classic(base_size = 10)+ ylab("calc_rmsd") +
  ggtitle("all frames") +
  xlab("ce_stage")



grid.arrange(p1,p2,p3,p4, nrow=2, ncol=2)
grid.arrange(p5,p6,p7,p8, nrow=2, ncol=2)
grid.arrange(p9,p10,p11,p12, nrow=2, ncol=2)

dev.off()


#cor.test(data.ts.sort$calc_rmsd, data.ts.sort$ce_score, method = "pearson")

## Output the pdb files corresponding to top X percent ce scores.
ext.dcd = paste(cohort,"*.dcd", sep = "/")

dcdfile = Sys.glob(ext.dcd)
# pdbfile = "/Users/bhandariyr/Documents/projects/afm/cobalamin/miscen/scripts/stats/contact10_5/pdb/cgRNaseP_FL_ori_center.pdb"
# pdbfile = args[1]
# dcdfile = args[2]
dcd <- read.dcd(dcdfile)
pdb <- read.pdb(pdbfile)
# ca.inds <- atom.select(pdb)
#setwd = "/Users/bhandariyr/Documents/projects/afm/cobalamin/miscen/scripts/stats/contact10_5"

for (val in data.ts.top$frame){
  if (val != 0) {
    pdb.name = paste("best_cc_frame_",val, ".pdb", sep = "")
    pdb.path = paste(dir_name, pdb.name, sep = "/")
    write.pdb(pdb, dcd[val, ], file=pdb.path)
  }
}

file1 = paste(dir_name, "energy_all.txt", sep = "/")
file2 = paste(dir_name, "/energy_top_", top_percent, "percent.txt", sep = "")
write.table(format(data.ts.sort), file = file1, sep = "\t", row.names = FALSE, quote = FALSE)
write.table(format(data.ts.top), file = file2, sep = "\t", row.names = FALSE, quote = FALSE)

print(paste("Correlation coef. of rmsd vs ce_score for all frames ",
      cor(data.ts.sort$calc_rmsd, data.ts.sort$ce_score, method = "pearson"), sep = ""))
print(paste("Correlation coef. of rmsd vs ce_score for top_", top_percent, "_percent_frames: ",
      cor(data.ts.top$calc_rmsd, data.ts.top$ce_score, method = "pearson"), sep = ""))




