## Usage: Rscript --vanilla analysis_bio3d.R pdbfile dcdfile time_stepfile atom_type
## atom_type = choose one (sugar, phosphate, base)
## the output will be text files afmcc.txt, rmsd.txt, rmsf.txt and plot file bio3d_plots.pdf

#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
#print(args[1])
#length(args)
## test if there is at least one argument: if not, return an error
if (length(args)!=4) {
  stop("Please provide the pdb, dcd, timestep file atom_type", call.=FALSE)
}
  
library(dplyr)
library(stringr)
library(bio3d)

# file = "/Users/bhandariyr/Documents/projects/afm/cafemol/dimer/twochain/run2_rescaled/output1/cbl_afm.ts"
file = args[3]
data.ts.lines = readLines(file)[grep("#all", readLines(file))]
data.ts =  data.ts.lines %>%
              strsplit("\\s+") %>%
              do.call(rbind, .) %>%
              data.frame()
dim(data.ts)
#View(data.ts)
data.ts$frame = as.numeric(rownames(data.ts)) - 1
data.ts.sort = data.ts[,c("frame", "X2", "X15")]
colnames(data.ts.sort) = c("frame","time.step", "afmcc")

## Name the pdf file to plot
out.file = "bio3d_plots.pdf"
pdf(out.file) 
plot(data.ts.sort[2:nrow(data.ts.sort),]$frame, data.ts.sort[2:nrow(data.ts.sort),]$afmcc, xlab = "Frame No.", ylab = "correlation")
write.table(format(data.ts.sort), file = "afmcc.txt", sep = "\t", row.names = FALSE, quote = FALSE)

# 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[4])
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
#rd <- rmsd(xyz[1,ca.inds$xyz], xyz[,ca.inds$xyz])
#rmsd.frame = data.frame(frame = 1:length(rd), rmsd = rd)
#write.table(rmsd.frame, file = "rmsd.txt", sep = "\t", row.names = FALSE, quote = FALSE)
#plot(rd, typ="l", ylab="RMSD", xlab="Frame No.")
#points(lowess(rd), typ="l", col="red", lty=2, lwd=2)
#hist(rd, breaks=40, freq=FALSE, main="RMSD Histogram", xlab="RMSD", ylab = "population")
#lines(density(rd), col="gray", lwd=3)
## root mean square flucturation
rf <- rmsf(xyz[,ca.inds$xyz])
rmsf.frame = data.frame(residue = 1:length(rf), rmsf = rf)
write.table(rmsf.frame, file = "rmsf.txt", sep = "\t", row.names = FALSE, quote = FALSE)
plot(rf, ylab="RMSF", xlab="Residue Position", typ="l")
## PCA analysis
# pc <- pca.xyz(xyz[,ca.inds$xyz])
# plot(pc, col=bwr.colors(nrow(xyz)), use.svd = TRUE )
## Cross correlation analysis
cij<-dccm(xyz[,ca.inds$xyz])
plot(cij)

# Close the pdf file
dev.off() 
