## Usage: Rscript --vanilla analysis_bio3d.R pdbfile dcdfile atom_type
## atom_type = choose one (sugar, phosphate, base, all)
## the output will be text files rmsd.txt

#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
if (length(args) != 3) {
  cat("Find RMSDs of all PDBs in the dcd file compared to any given PDB\n")
  cat("Usage:\n")
  cat("       Rscript get_rmsd.R pdb_file dcd_file atom_type\n")
  cat("       pdb_file     reference PDB file compare against\n")
  cat("       dcd_file     trajectory file\n")
  cat("       atom_type    choose one (sugar, phosphate, base, all)\n")
  quit()
}

  
library(dplyr)
library(stringr)
library(bio3d)

# dcdfile = "~/Documents/projects/afm/cafemol/dimer/twochain/run2_rescaled/output1/cbl_afm.dcd"
# pdbfile = "~/Documents/projects/afm/cafemol/dimer/twochain/run1/pdb/cbl_dimer_offset.pdb"
pdbfile = args[1]
dcdfile = args[2]
dcd <- read.dcd(dcdfile)
pdb <- read.pdb(pdbfile)

atom.type = as.character(args[3])
if (atom.type == "sugar") {
   ca.inds <- atom.select(pdb, elety = "S")
}   else if (atom.type == "phosphate") {
    ca.inds <- atom.select(pdb, elety = "P")
}  else if (atom.type == "base") {
    ca.inds <- atom.select(pdb, elety = c("Ab","Ub", "Gb","Cb"))
}  else if (atom.type == "all") {
    ca.inds <- atom.select(pdb)
}

ca.inds
# ca.inds <- atom.select(pdb)
# Print details of the first few selected atoms
# pdb$atom[ca.inds$atom, ]

xyz <- fit.xyz(fixed=pdb$xyz, mobile=dcd,
               fixed.inds=ca.inds$xyz,
               mobile.inds=ca.inds$xyz)
# dim(xyz) == dim(dcd)

## Root mean square deviation
print("Now calculating rmsd")
#rd <- rmsd(xyz[1,ca.inds$xyz], xyz[,ca.inds$xyz])
rd <- rmsd(pdb$xyz[1, ca.inds$xyz], xyz[,ca.inds$xyz])
rmsd.frame = data.frame(frame = 1:length(rd), rmsd = rd)
out.file = paste0("rmsd_", atom.type,".txt",sep="")
write.table(rmsd.frame, file = out.file, sep = "\t", row.names = FALSE, quote = FALSE)
