#' Plots vaf distribution of genes
#' @description Plots vaf distribution of genes as a boxplot.
#'
#' @param maf an \code{\link{MAF}} object generated by \code{\link{read.maf}}
#' @param vafCol manually specify column name for vafs. Default looks for column 't_vaf'
#' @param genes specify genes for which plots has to be generated
#' @param orderByMedian Orders genes by decreasing median VAF. Default TRUE
#' @param top if \code{genes} is NULL plots top n number of genes. Defaults to 5.
#' @param flip if TRUE, flips axes. Default FALSE
#' @param showN if TRUE, includes number of observations
#' @param gene_fs font size for gene names. Default 0.8
#' @param fn Filename. If given saves plot as a output pdf. Default NULL.
#' @param axis_fs font size for axis. Default 0.8
#' @param width Width of plot to be saved. Default 4
#' @param height Height of plot to be saved. Default 5
#' @return ggplot object which can be further modified.
#' @examples
#' laml.maf <- system.file("extdata", "tcga_laml.maf.gz", package = "maftools")
#' laml <- read.maf(maf = laml.maf)
#' plotVaf(maf = laml, vafCol = 'i_TumorVAF_WU')
#'
#' @export

plotVaf = function(maf, vafCol = NULL, genes = NULL, top = 10,
                   orderByMedian = TRUE, flip = FALSE, fn = NULL,
                   gene_fs = 0.8, axis_fs = 0.8, height = 5, width = 5, showN = TRUE){

  if(is.null(genes)){
    genes = as.character(getGeneSummary(x =maf)[1:top, Hugo_Symbol])
  }

  dat = subsetMaf(maf = maf, genes = genes, includeSyn = FALSE)

  if(!'t_vaf' %in% colnames(dat)){
    if(is.null(vafCol)){
      if(all(c('t_ref_count', 't_alt_count') %in% colnames(dat))){
        message("t_vaf field is missing, but found t_ref_count & t_alt_count columns. Estimating vaf..")
        dat[,t_vaf := as.numeric(as.character(t_alt_count))/(as.numeric(as.character(t_ref_count)) + as.numeric(as.character(t_alt_count)))]
      }else{
        print(colnames(dat))
        stop('t_vaf field is missing. Use vafCol to manually specify vaf column name.')
      }
    }else{
      colnames(dat)[which(colnames(dat) == vafCol)] = 't_vaf'
    }
  }

  #dat.genes = data.frame(dat[dat$Hugo_Symbol %in% genes])
  #suppressMessages(datm <- melt(dat.genes[,c('Hugo_Symbol', 't_vaf')]))
  dat.genes = dat[dat$Hugo_Symbol %in% genes]
  suppressWarnings(datm <- data.table::melt(data = dat.genes[,.(Hugo_Symbol, t_vaf)]))
  #remove NA from vcf
  datm = datm[!is.na(value)]

  #maximum vaf
  if(max(datm$value, na.rm = TRUE) > 1){
    datm$value = datm$value/100
  }

  if(orderByMedian){
    geneOrder = datm[,median(value),Hugo_Symbol][order(V1, decreasing = TRUE)][,Hugo_Symbol]
    datm$Hugo_Symbol = factor(x = datm$Hugo_Symbol, levels = geneOrder)
  }


  bcol = c(RColorBrewer::brewer.pal(n = 8, name = "Dark2"),
           RColorBrewer::brewer.pal(n = 8, name = "Accent"))
  if(length(genes) > length(bcol)){
    bcol = sample(x = colors(distinct = TRUE), size = length(genes), replace = FALSE)
  }

  if(!is.null(fn)){
    pdf(file = paste0(fn, ".pdf"), width = width, height = height, paper = "special", bg = "white")
  }

  if(flip){
    par(mar = c(3, 4, 2, 2))
    b = boxplot(value ~ Hugo_Symbol, data = datm, xaxt="n", boxwex=0.5, outline=FALSE, lty=1, lwd = 1.4, outwex=0,
                staplewex=0, ylim = c(0, 1), axes = FALSE, border = bcol, horizontal = TRUE)
    axis(side = 1, at = seq(0, 1, 0.2), las = 1, font = 2, lwd = 2.2, cex.axis = axis_fs)
    axis(side = 2, at = 1:length(b$names), labels = b$names, tick = FALSE, las = 2, font = 2, line = -1, cex.axis = gene_fs)
    if(showN){
      axis(side = 4, at = 1:length(b$names), labels = b$n, font = 2, tick = FALSE, line = -1, las = 2, cex.axis = gene_fs)
    }


    stripchart(value ~ Hugo_Symbol, vertical = FALSE, data = datm,
               method = "jitter", add = TRUE, pch = 16, col = bcol, cex = 0.5)

  }else{
    par(mar = c(4, 3, 2, 1))
    b = boxplot(value ~ Hugo_Symbol, data = datm, xaxt="n", boxwex=0.5, outline=FALSE, lty=1, lwd = 1.4, outwex=0,
                staplewex=0, ylim = c(0, 1), axes = FALSE, border = bcol)
    axis(side = 1, at = 1:length(b$names), labels = b$names, tick = FALSE, las = 2, font = 2, line = -1, cex.axis = gene_fs)
    axis(side = 2, at = seq(0, 1, 0.2), las = 2, cex.axis = axis_fs, lwd = 3, font.axis = 2, cex = 1.5, font = 2)
    if(showN){
      axis(side = 3, at = 1:length(b$names), labels = b$n, font = 2, tick = FALSE,
           line = -1, cex.axis = gene_fs)
    }

    stripchart(value ~ Hugo_Symbol, vertical = TRUE, data = datm,
               method = "jitter", add = TRUE, pch = 16, col = bcol, cex = 0.5)

  }

  if(!is.null(fn)){
    dev.off()
  }
}
