#' Plot differences between APOBEC enriched and non-APOBEC enriched samples.
#'
#' @description Plots differences between APOBEC enriched and non-APOBEC enriched samples
#'
#' @details Plots differences between APOBEC enriched and non-APOBEC enriched samples (TCW). Plot includes differences in mutations load, tCw motif distribution and top genes altered.
#'
#' @param tnm output generated by \code{\link{trinucleotideMatrix}}
#' @param maf an \code{\link{MAF}} object used to generate the matrix
#' @param title_size size of title. Default 1.3
#' @param axis_lwd axis width. Default 2.
#' @param font_size font size. Default 1.2
#' @return list of table containing differenatially altered genes. This can be passed to \code{\link{forestPlot}} to plot results.
#' @examples
#' \dontrun{
#' laml.tnm <- trinucleotideMatrix(maf = laml, ref_genome = 'hg19.fa', prefix = 'chr',
#' add = TRUE, useSyn = TRUE)
#' plotApobecDiff(laml.tnm)
#' }
#' @importFrom grDevices pdf boxplot.stats dev.off
#' @seealso \code{\link{trinucleotideMatrix}} \code{\link{plotSignatures}}
#' @export

plotApobecDiff = function(tnm, maf, title_size = 1, axis_lwd = 2, font_size = 1.2){
  sub.tbl <- tnm$APOBEC_scores
  sub.tbl$APOBEC_Enriched = factor(sub.tbl$APOBEC_Enriched, levels = c('yes', 'no')) #Set levels
  yp = c(min(boxplot.stats(x = sub.tbl[,n_mutations])$stats), max(boxplot.stats(x = sub.tbl[,n_mutations])$stats))

  if(length(boxplot.stats(x = sub.tbl[,n_mutations])$out) > 0){
    yp[2] = min(boxplot.stats(x = sub.tbl[,n_mutations])$out)+3
  }else{
    yp[2] = yp[2]+3
  }

  yp = as.integer(seq(yp[1], yp[2], length.out = 5))
  if(yp[1] <= 3){
    yp[1] = 0
  }else{
    yp[1] = yp[1] - 3
  }

  if(nrow(sub.tbl[!is.na(APOBEC_Enriched), mean(fraction_APOBEC_mutations), APOBEC_Enriched][APOBEC_Enriched %in% 'yes']) == 0){
    stop('None of the samples are enriched for APOBEC. Nothing to plot.')
  }

  apobec.samps = as.character(sub.tbl[APOBEC_Enriched %in% 'yes', Tumor_Sample_Barcode])
  apobec.maf = subsetMaf(maf = maf, tsb = apobec.samps, mafObj = TRUE)
  apobec.maf.summary = apobec.maf@summary[ID %in% 'total', .(Mean, Median)][,Cohort := 'Enriched']

  nonapobec.samps = as.character(sub.tbl[APOBEC_Enriched %in% 'no', Tumor_Sample_Barcode])
  nonapobec.maf = subsetMaf(maf = maf, tsb = nonapobec.samps, mafObj = TRUE)
  nonapobec.maf.summary = nonapobec.maf@summary[ID %in% 'total', .(Mean, Median)][,Cohort := 'nonEnriched']

  mc = mafCompare(m1 = apobec.maf, m2 = nonapobec.maf, m1Name = "Enriched", m2Name = "nonEnriched", minMut = 2, useCNV = TRUE)
  mcr = mc$results
  mc$SampleSummary = merge(mc$SampleSummary, rbind(apobec.maf.summary, nonapobec.maf.summary))

  if(nrow(mcr[pval < 0.05]) < 10){
    mcg = mcr[1:10, Hugo_Symbol]
  }else{
    mcg = mcr[pval < 0.05, Hugo_Symbol][1:10]
  }

  apobec.mcg.gs = getGeneSummary(x = apobec.maf)[Hugo_Symbol %in% mcg, .(Hugo_Symbol, MutatedSamples)]
  colnames(apobec.mcg.gs)[2] = 'Enriched'
  nonapobec.mcg.gs = getGeneSummary(x = nonapobec.maf)[Hugo_Symbol %in% mcg, .(Hugo_Symbol, MutatedSamples)]
  colnames(nonapobec.mcg.gs)[2] = 'nonEnriched'

  mcg.gs = merge(apobec.mcg.gs, nonapobec.mcg.gs, all = TRUE)
  mcg.gs[is.na(mcg.gs)] = 0

  mcg.gs$Enriched = round(mcg.gs$Enriched / mc$SampleSummary[Cohort %in% colnames(mcg.gs)[2], SampleSize], digits = 3)
  mcg.gs$nonEnriched = round(mcg.gs$nonEnriched / mc$SampleSummary[Cohort %in% colnames(mcg.gs)[3], SampleSize], digits = 3)

  data.table::setDF(mcg.gs, rownames = mcg.gs$Hugo_Symbol)
  mcg.gs = as.matrix(mcg.gs[,-1])
  mcg.gs = mcg.gs[mcg,, drop = FALSE]

  pieDat = sub.tbl[!is.na(APOBEC_Enriched), mean(fraction_APOBEC_mutations), APOBEC_Enriched]
  pieDat[,nonApobec := 1 - V1]
  colnames(pieDat)[2] = 'Apobec'
  pieDat = data.table::melt(pieDat, id.vars = 'APOBEC_Enriched', drop = FALSE)
  pieDat[,title := paste0(variable, ' [', round(value, digits = 3), ']')]
  #pieDat$title = gsub(pattern = '^Apobec', replacement = 'tCw', x = pieDat$title)
  #pieDat$title = gsub(pattern = '^nonApobec', replacement = 'non-tCw', x = pieDat$title)
  pieDat$title = round(pieDat$value, digits = 3)

  pieCol  = c("#084594", "#9ECAE1")

  layout(matrix(c(1,2,3,1,4,4), 2, 3, byrow = TRUE), widths=c(2, 2, 2))
  par(bty="n", mgp = c(0.5,0.5,0), mar = c(2,3.5,4,0) + 0.1, las=1, tcl=-.25, font.main=4, xpd=NA)

  boxplot(at = 1:2, n_mutations ~ APOBEC_Enriched, data = sub.tbl,  xaxt="n", boxwex=0.6, outline = TRUE, lty=1,
          outwex=0, staplewex=0, frame.plot = FALSE, col = c('maroon', 'royalblue'), yaxt = 'n',
          ylim = c(min(yp), max(yp)),
          outcol="gray70", outcex = 0.8, outpch  = 16, boxfill = NULL,
          border = c('maroon', 'royalblue'), lwd = 1.6)

  title(main = 'Mutation load between APOBEC &\nnon-APOBEC enriched samples', cex.main = title_size)

  axis(side = 1, at = c(1, 2), labels = na.omit(sub.tbl[,.N,APOBEC_Enriched])[,paste0('N=', N)],
       las = 1, tick = FALSE, font.axis = 2, cex = 2.5, font = 2, cex.axis = font_size, lwd = axis_lwd)
  axis(side = 2, at = yp, lwd = axis_lwd, las = 1, font = 2, cex.axis = font_size, font.axis = 2)

  p = wilcox.test(sub.tbl[APOBEC_Enriched %in% 'yes', n_mutations], sub.tbl[APOBEC_Enriched %in% 'no', n_mutations])$p.value
  if(p < 0.001 ){
    #lines(x = c(1.3, 1.8), y = rep(yp[length(yp)]-3, 2), lwd = 2)
    text(x = 1.55, y = yp[length(yp)], labels = "***", cex = 3)
  }else if(p < 0.01){
    text(x = 1.55, y = yp[length(yp)], labels = "**", cex = 3)
  }else if(p < 0.05){
    text(x = 1.55, y = yp[length(yp)], labels = "*", cex = 3)
  }


  par(bty="n", mgp = c(0.5,0.5,0), mar = c(0.5,1.5,3,1) + 0.1, las=1, tcl=-.25, font.main=4, xpd=NA)
  pie(x = pieDat[APOBEC_Enriched %in% 'yes', value], col = pieCol,
      border="white", radius = 0.95, cex.main=0.6, cex = font_size, font = 2,
      labels =  round(pieDat[APOBEC_Enriched %in% 'yes', title], digits = 2), clockwise = TRUE)
  symbols(0,0,circles=.4, inches=FALSE, col="white", bg="white", lty=0, add=TRUE)
  title(main = 'tCw load\nAPOBEC samples', cex.main = title_size)

  pie(x = pieDat[APOBEC_Enriched %in% 'no', value], col = pieCol,
      border="white",  radius = 0.95, cex.main = 1.33, cex = font_size, font = 2,
      labels =  round(pieDat[APOBEC_Enriched %in% 'no', title], 2), clockwise = TRUE)
  symbols(0,0,circles=.4, inches=FALSE, col="white", bg="white", lty=0, add=TRUE)
  title(main = 'tCw load\nnon-APOBEC samples', cex.main = title_size)

  rownames(mcg.gs) = paste0(rownames(mcg.gs), ifelse(test = mcr[Hugo_Symbol %in% rownames(mcg.gs), pval] < 0.001, yes = '***',
                                  no = ifelse(test = mcr[Hugo_Symbol %in% rownames(mcg.gs), pval] < 0.01, yes = '**',
                                              no = ifelse(mcr[Hugo_Symbol %in% rownames(mcg.gs), pval] < 0.05, yes = '*', no = ''))))

  par(bty="n", mgp = c(0.5,0.5,0), mar = c(7,3.5,0,0) + 0.1, las=1, tcl=-.25, font.main=4, xpd=NA)
  barplot(height = t(mcg.gs), beside = TRUE, las =2, ylab = "", border = NA, font.axis = 2,
          font.lab = 2, col = c('maroon', 'royalblue'),
          ylim = c(0, yl = max(apply(mcg.gs, 2, max))), yaxt = "n", cex.names = font_size)
  axis(side = 2, at = seq(0, max(apply(mcg.gs, 2, max)), length.out = 4), lwd = axis_lwd, font.axis = 2, cex = 1.5,
        labels = paste0(as.integer((seq(0, max(apply(mcg.gs, 2, max)), length.out = 4))*100, digits = 2), '%'), cex.axis = font_size)

  return(mc)
}
