#!/usr/bin/Rscript

suppressMessages(library('plotly'))
suppressMessages(library('heatmaply'))

#packageVersion('plotly')

get_args <- function(){
  
  ## Collect arguments
  args <- commandArgs(TRUE)
  
  ## Default setting when no arguments passed
  if(length(args) < 1) {
    args <- c("--help")
  }
  
  ## Help section
  if("--help" %in% args) {
    cat("Pathview R script
    Arguments:
      --help                  Print this test
      --input                 path of the input  file (must contains a colum of uniprot and/or geneID accession number)
      --output                Output file
      --type                  type of output file, could be html, pdf, jpg or png
      --cols                  Columns to use for heatmap, exemple : '3:8' to use columns from the third to the 8th
      --row_names             Column which contains row names
      --header                True or False
      --col_text_angle        Angle of columns label ; from -90 to 90 degres

      Example:
      ./heatmap_viz.R --input='dat.nucl.norm.imputed.tsv' --output='heatmap.html' --cols='3:8' --row_names='2' --header=TRUE --col_text_angle=0 \n\n")
    
    q(save="no")
  }
  
  parseArgs <- function(x) strsplit(sub("^--", "", x), "=")
  argsDF <- as.data.frame(do.call("rbind", parseArgs(args)))
  args <- as.list(as.character(argsDF$V2))
  names(args) <- argsDF$V1
  
  return(args)
}

read_file <- function(path,header){
  file <- try(read.table(path,header=header, sep="\t",stringsAsFactors = FALSE, quote="",fill=TRUE),silent=TRUE)
  if (inherits(file,"try-error")){
    stop("File not found !")
  }else{
    return(file)
  }
}

#convert a string to boolean
str2bool <- function(x){
  if (any(is.element(c("t","true"),tolower(x)))){
    return (TRUE)
  }else if (any(is.element(c("f","false"),tolower(x)))){
    return (FALSE)
  }else{
    return(NULL)
  }
}

#remove remaining quote 
#remove lines with at least one empty cell in a matrix between two defined columns
clean_df <- function(mat,first_col,last_col,rownames){
  tmp = mat[,first_col:last_col]
  tmp <- as.data.frame(apply(tmp,c(1,2),function(x) {ifelse(is.character(x),as.numeric(x),x)}))
  bad_lines <- which(apply(tmp, 1, function(x) any(is.na(x))))
  mat <- cbind(mat[,as.numeric(rownames)],tmp)
  if (length(bad_lines) > 0) {
    mat <- mat[- bad_lines,]
    print(paste("lines",bad_lines, "has been removed: at least one non numeric content"))
  }
  return(mat)
}

#get args
args <- get_args()

header=str2bool(args$header)
output <- rapply(strsplit(args$output,"\\."),c) #remove extension
output <- paste(output[1:length(output)-1],collapse=".")
output <- paste(output,args$type,sep=".")
first_col=as.numeric(substr(args$cols,1,1))
last_col=as.numeric(substr(args$cols,3,3))

#cleaning data
uto <- read_file(args$input,header = header)
uto <- clean_df(uto,first_col,last_col,args$row_names)
data <- as.data.frame(uto[,-1])
row_names = uto[,1]
if (header) {
  col_names = names(data)
} else {
  col_names = c(first_col:last_col)
}

#building heatmap
heatmaply(data, file=output, margins=c(100,50,NA,0), plot_method="plotly", labRow = row_names, labCol = col_names,
          grid_gap = 0,cexCol = 1, column_text_angle = as.numeric(args$col_text_angle), width = 1000, height=1000, colors = c('blue','green','yellow','red'))


####heatmaply

simulateExprData <- function(n, n0, p, rho0, rho1){ 
  # n: total number of subjects 
  # n0: number of subjects with exposure 0 
  # n1: number of subjects with exposure 1 
  # p: number of genes 
  # rho0: rho between Z_i and Z_j for subjects with exposure 0 
  # rho1: rho between Z_i and Z_j for subjects with exposure 1
  
  # Simulate gene expression values according to exposure 0 or 1, 
  # according to a centered multivariate normal distribution with 
  # covariance between Z_i and Z_j being rho^|i-j| 
  n1 <- n - n0 
  times <- 1:p
  H <- abs(outer(times, times, "-")) 
  V0 <- rho0^H 
  V1 <- rho1^H 
  
  # rows are people, columns are genes 
  genes0 <- MASS::mvrnorm(n = n0, mu = rep(0,p), Sigma = V0) 
  genes1 <- MASS::mvrnorm(n = n1, mu = rep(0,p), Sigma = V1) 
  genes <- rbind(genes0,genes1) 
  return(genes)
}

#genes <- simulateExprData(n = 50, n0 = 25, p = 100, rho0 = 0.01, rho1 = 0.95)

#heatmaply(genes, k_row = 2, k_col = 2)

#heatmaply(cor(genes), k_row = 2, k_col = 2)