## Process the image data produced by cafemol
## Usage: Rscript cafemol_image.txt cafemol_input.txt pdb.name

library(dplyr)
library(stringr)
#!/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)!=3) {
  stop("Please provide cafemol_image.txt cafemol_input.txt pdb.name", call.=FALSE)
}


# image.file = "/Users/bhandariyr/Documents/projects/afm/cafemol/miscen/scripts/image_info_cafe.txt"
image.file = args[1]
df = read.table(image.file, 
               sep="",
               fill=FALSE, 
               strip.white=TRUE)
# dim(df) 
# head(df)
colnames(df)= c("serial", "rank", "expt", "calc") 
df1 = df[1:(nrow(df)/2), ]
# dim(df1)
df1.sort = df1[order(df1$rank), ]
# View(df1.sort)
df2 = df[(nrow(df)/2+1):nrow(df), ]
df2.sort = df2[order(df2$rank), ]
# View(df2.sort)
## Calculate the correlation coefficient
cc1 = sum(df1$expt*df1$calc)/(sqrt(sum(df1$expt**2))*sqrt(sum(df1$calc**2)))
cc2 = sum(df2$expt*df2$calc)/(sqrt(sum(df2$expt**2))*sqrt(sum(df2$calc**2)))
# cc1; cc2

## read cafemol input file to validate the cafemol results
## and to get the template for x and y pixels
# input.file = "/Users/bhandariyr/Documents/projects/afm/cafemol/miscen/scripts/cbl_afm_fit1.inp"
input.file = args[2]
data.input.lines = readLines(input.file)[grep("AFM_IMAGE", readLines(input.file))]
data.cafeinput =  data.input.lines %>%
  strsplit("\\s+") %>%
  do.call(rbind, .) %>%
  data.frame()
# head(data.cafeinput)
# class(data.cafeinput); dim(data.cafeinput)
colnames(data.cafeinput) = c("name","x.pos", "y.pos", "z")
# plot(df1.sort$expt, data.cafeinput$z)
# head(df1.sort$expt)
# head(data.cafeinput$z)

# out.file = "/Users/bhandariyr/Documents/projects/afm/cafemol/miscen/scripts/out.txt"
out.file = paste(args[3], "calc_data.txt", sep = "_")
corr.info = paste("AFMcc of pdb file ", args[3], " is ", cc1 )
corr.info
write.table(cbind(data.cafeinput$x.pos, data.cafeinput$y.pos, df1.sort$calc), out.file,
            append = FALSE, sep = " ", row.names = FALSE, col.names = c("x","y","z"), quote = FALSE)






