# Helper functions
# fdrcut function creates a new factor column based on specified numeric cuts on a p-value column for plotting
fdrcut <- function(adjpval, fdr.range.cut =  c(0,0.02,0.05,1)) {
  
  adjpval.ranges <- cut(adjpval, fdr.range.cut, include.lowest = T)  
  fdr.range.cut <- sprintf("%.2f", fdr.range.cut)
  if(length(fdr.range.cut) == 4) {
    adjpval.ranges <- factor(adjpval.ranges,levels = rev(levels(adjpval.ranges)),ordered = T,
                             labels = rev( c( paste(fdr.range.cut[1],'< q <=',fdr.range.cut[2],sep = ''),
                                              paste(fdr.range.cut[2],'< q <=',fdr.range.cut[3],sep = ''),
                                              paste(fdr.range.cut[3],'< q <=',fdr.range.cut[4],sep = ''))))
  } else if(length(fdr.range.cut) == 5) {
    adjpval.ranges <- factor(adjpval.ranges,levels = rev(levels(adjpval.ranges)),ordered = T,
                             labels = rev( c( paste(fdr.range.cut[1],'< q <=',fdr.range.cut[2],sep = ''),
                                              paste(fdr.range.cut[2],'< q <=',fdr.range.cut[3],sep = ''),
                                              paste(fdr.range.cut[3],'< q <=',fdr.range.cut[4],sep = ''),
                                              paste(fdr.range.cut[4],'< q <=',fdr.range.cut[5],sep = ''))))
    
  } else {
    stop("fdr.range.cut vector not of length 4")
  }
  return(adjpval.ranges)
}

# return a table of differentially expressed proteins from an eBayes fitted linear model object
get.toptable <- function(fitobj, coef, glist, suffix) {
  ttemp <- topTable(fitobj,coef = coef,genelist = glist,number = Inf,p.value = 1)
  datinds <-  which(names(ttemp) != "ID")
  names(ttemp)[datinds] <- paste(names(ttemp)[datinds], suffix, sep="")
  return(ttemp) 
}

# Helper function to create an interactive MA plot using the R Plotly package. Make sure gene symbols in tdf are indicated in column Gene.name
ggmaplot <- function(tdf, # data frame containing the data to plot
                 xvar, # string name of x-variable
                 yvar,# string name of y-variable
                 adjpvalvar,# string name of adjusted p-value to color points by significance
                 fdr.range.cut = c(0,0.01,0.05,0.1,1), # cuts used to categorize adj.p.vals
                 high.sig.threshold = 0.0000001, # all points with adj.p.value <= to this are annotated with gene symbols that must be named Gene.name, 
                 xlimits = Inf,
                 ylimits = Inf,
                 xtitle = NULL,
                 ytitle = NULL,
                 plottitle = NULL,
                 annotgenes = T,
                 revcols = T) {
  require(ggplot2)
  require(RColorBrewer)
  require(plotly)
  source('D:/Will/R/Functions/RSparta/R/fdrfunctions.R')
#  browser()
  xvar <- enquo(xvar)
  yvar <- enquo(yvar)
  adjpvalvar <- enquo(adjpvalvar) 
  adjpvalvarcut <- fdrcut(tdf %>% select(!! adjpvalvar) %>% pull(), fdr.range.cut = fdr.range.cut)
  
  if(revcols) {
    rcols <- rev(brewer.pal(n = 7, name = "Reds"))[1:(length(fdr.range.cut) - 1)]
  } else {
    rcols <- brewer.pal(n = 7, name = "Reds")[1:(length(fdr.range.cut) - 1)]
  }
  
  gcol <- brewer.pal(n = 5, name = "Greys")[3]
  
  ggy <- tdf %>% ggplot(mapping = aes_string( quo_name(xvar), quo_name(yvar))) +
    geom_point(aes(text = Gene.name,
                   color = adjpvalvarcut)) +
    scale_colour_manual(values = c(gcol, rev(rcols))) +
    geom_vline(xintercept = 0,color = "grey") +
    theme_classic(base_size = 15,
                  base_line_size = 2)  +
     theme(legend.title = element_blank())
 
  if ( !is.infinite(xlimits[1])) {
    ggy <- ggy + xlim(xlimits)
  } 
  if ( !is.infinite(ylimits[1])) {
    ggy <- ggy + ylim(ylimits)
  } 
  if(!is.null(xtitle) ) {
    ggy <- ggy + xlab(xtitle)
  }
  
  if(!is.null(ytitle) ) {
    ggy <- ggy + xlab(ytitle)
  }
  
  if(!is.null(plottitle) ) {
    ggy <- ggy +  ggtitle(plottitle) 
  }
  
  tdf.thresh10pct <- tdf %>% filter(!! adjpvalvar <= 0.01 & !! adjpvalvar > 0.0001)
    
  tdf.thresh1pct <- tdf %>% filter(!! adjpvalvar <= high.sig.threshold)
  
  tdf.threshy <- tdf %>% filter(!! yvar > 25 & !! xvar < 2 )  

  ggply <- ggplotly(ggy, tooltip = "text")
  
  ggply <- ggply %>% add_annotations(x = tdf.thresh10pct %>% select(!! xvar) %>% pull(),
                                     y  = tdf.thresh10pct %>% select(!! yvar) %>% pull(),
                                     text = tdf.thresh10pct$Gene.name,
                                     font = list(family = 'Arial',
                                                 size = 12,
                                                 color = 'black'),
                                     showarrow = T,
                                     visible = F, 
                                     arrowhead = 4,
                                     arrowwidth = 0.5,
                                     arrowcolor = 'grey',
                                     arrowsize = .5,
                                     ax = 20,
                                     ay = -20,
                                     clicktoshow = "onoff")
  if (annotgenes) {
    ggply <- ggply %>% add_annotations(x = tdf.thresh1pct %>% select(!! xvar) %>% pull(),
                    y  = tdf.thresh1pct %>% select(!! yvar) %>% pull(),
                    text = tdf.thresh1pct$Gene.name,
                    font = list(family = 'Arial',
                                size = 10,
                                color = 'black'),
                    showarrow = T,
                    visible = T, 
                    arrowhead = 4,
                    arrowwidth = 0.5,
                    arrowcolor = 'grey',
                    arrowsize = .5,
                    ax = 20,
                    ay = -20,
                    clicktoshow = "onoff") 
  }
  
  ggply <- ggply %>% add_annotations(x = tdf.threshy %>% select(!! xvar) %>% pull(),
                  y  = tdf.threshy %>% select(!! yvar) %>% pull(),
                  text = tdf.threshy$Gene.name,
                  font = list(family = 'Arial',
                              size = 12,
                              color = 'black'),
                  showarrow = T,
                  visible = F, 
                  arrowhead = 4,
                  arrowwidth = 0.5,
                  arrowcolor = 'grey',
                  arrowsize = .5,
                  ax = 20,
                  ay = -20,
                  clicktoshow = "onoff")
  
 
  
  
  return(list(maplot.ggplot=ggy,
              maplot.plotly = ggply))

  
}

# Try:
# myggma <- tophits %>% ggmaplot(xvar = abcam, yvar = avgabcam, adjpvalvar = abcam.adj.P.Val)
# ggplotly(myggma, tooltip = "text")
