dauergut: Scripts for rict-1 dauer paper

#' plot_CIs
#'
#' Function generates boxplots for dauer, expression or roaming data along with
#' accompanying NHST comparisons and 95 pct confidence intervals
#' @param df input dataset. Requires a "genotype" column. Response needs to be "pct" for
#' dauer plots, "n_entries" for roaming data, and "cell.norm" for gene expression microscopy
#' data. Also requires a "mixed" dataset, generated by either MM_CIs.R, Mm_CI_trans_Tukey.R or MM_CI_trans_dunnett.R
#' @param title title for the plot
#' @param plot.contrasts list of pairwise comparisons to plot. Usually generated by selecting a list 
#' of comparisons using tukey_contrasts or dunnett_contrasts functions which use the lsmeans
#' package, ie plot.contrasts = c("",contrasts$prange[1:2]) for 2
#' comparisons to a control.
#' @param plot.contrasts.2 list of secondary comparisons (ie rescue lines). These will be colored red and offset in the
#' y direction indicated by offset.
#' @param ypos y position of p value for plot.contrasts on the plot response scale
#' @param type type of plot ("dauer", "roam", "GFP"). Defaults to mRNA plot if no input.
#' @param offset position to offset plot.contrast.2 on the y response scale. 
#' @export
#' @examples plot_CIs(df, "GFP plot", plot.contrasts = contrasts, ypos = 500, type = "GFP", offset = 100)
plot_CIs<-function (df, title, plot.contrasts, plot.contrasts.2, ypos, type, offset, labels) {
  #generate plots for different types of data
  #all using genotype as predictor factor
  #type = "dauer" or "grid" or "GFP" or mRNA expression plot
  #set common values for all plots:
  box.width <- 0.3
  point.size <- 1
  line.width <- 0.15
  text.size <- 4
  p <- ggplot(df, aes(x=genotype)) + #x-layer
    theme_my +
    geom_errorbar(data=mixed, aes(x=x.pos,y=mean, ymin=lower.CL, ymax=upper.CL),
                  width=0,colour ="grey", lwd=line.width) + #90% cred int for 95% one-sided H0
    geom_errorbar(data=mixed, aes(x=x.pos,y=mean, ymin = lower.25, ymax = upper.75),
                  width=0,colour = "darkgrey", lwd = line.width+0.7) + #75% cred interval
    geom_segment(data = mixed, aes(x = x.pos-0.005*nrow(mixed),
                                   y = mean, xend = x.pos+0.005*nrow(mixed), 
                                   yend = mean), colour = "darkgrey") + 
    scale_x_discrete(labels=function(x) sub(" ","\n",x,fixed=TRUE)) +
    stat_summary(aes(x=genotype, y=ypos), geom="text", label=plot.contrasts, show.legend = TRUE, size=text.size) + # pvalues
    theme(axis.text.x = element_text(size = 12),
          axis.text.y = element_text(size = 12),
          axis.line = element_line(size=0.2),
          axis.title = element_text(size=16))
  if(type == "dauer") {
    p1 <- p + 
      geom_dotplot(aes(y=pct, x = genotype),binwidth=.015, binaxis="y", position="dodge", stackdir="center", size =.3) +
      stat_summary(aes(y=pct),fun.y = median, fun.ymin = median, fun.ymax = median,
                  geom = "crossbar", width = 0.25, lwd = line.width) +
      labs(title = title,
           subtitle="dauer plot",
           y = "proportion dauer",
           x = "genotype"
      ) +
      scale_y_continuous(breaks=c(0,0.25,0.5,0.75,1.0))
  } else {#roaming plots
    if(type == "grid") {
      p1 <- p + 
        #geom_boxplot(aes(y=n_entries), width=box.width, outlier.shape=NA, lwd=line.width, fill = "#33CCFF") +
        geom_dotplot(aes(y=n_entries), binaxis = "y",
                     binwidth = 7.5,
                     stackdir = "center",
                     colour = "#482677FF",
                     fill = "#482677FF") +
        # geom_quasirandom(aes(y=n_entries),colour = "grey", cex=1,
        #                  width = 0,size=0.3*point.size,
        #                  method = 'smiley') +
        #geom_point(aes(y=n_entries), size=0.3*point.size)) + 
        stat_summary(aes(y=n_entries),fun.y = median, 
                     fun.ymin = function(z) {quantile(z,0.25)}, 
                     fun.ymax = function(z) {quantile(z,0.75)},
                     geom = "errorbar", width = 0.15, lwd = line.width) +
        stat_summary(aes(y=n_entries),fun.y = median, 
                     fun.ymin = median, 
                     fun.ymax = median,
                     geom = "crossbar", width = 0.25, lwd = line.width+0.2) +
        #stat_boxplot(aes(y=n_entries), geom ='errorbar') +
        labs(title = title,subtitle="roaming plot",y = "grid entries",x = "genotype")
    } else {#GFP expression plots
      if(type == "GFP") {
        p1 <- p +
          geom_quasirandom(aes(y=cell.norm),colour = "#339900", cex=1,
                           width = 0.075,size=0.3*point.size,
                           method = 'smiley') +
          stat_summary(aes(y=cell.norm),fun.y = median, 
                       fun.ymin = function(z) {quantile(z,0.25)}, 
                       fun.ymax = function(z) {quantile(z,0.75)},
                       geom = "errorbar", width = 0.15, lwd = line.width) +
          stat_summary(aes(y=cell.norm),fun.y = median, 
                       fun.ymin = median, 
                       fun.ymax = median,
                       geom = "crossbar", width = 0.25, lwd = line.width+0.3) +
          labs(title = title,subtitle="GFP plot",y = "normalized expression",x = "genotype")
      } else {#mRNA plots
        p1 <- p +
          geom_quasirandom(aes(y=cell.norm),colour = "#990000", cex=1,
                           width = 0.075,size=0.3*point.size,
                           method = 'smiley') +
          stat_summary(aes(y=cell.norm),fun.y = median, 
                       fun.ymin = function(z) {quantile(z,0.25)}, 
                       fun.ymax = function(z) {quantile(z,0.75)},
                       geom = "errorbar", width = 0.15, lwd = line.width) +
          stat_summary(aes(y=cell.norm),fun.y = median, 
                       fun.ymin = median, 
                       fun.ymax = median,
                       geom = "crossbar", width = 0.25, lwd = line.width+0.2) +
          labs(title = title,subtitle="mRNA plot",y = "normalized expression",x = "genotype")
          # geom_boxplot(aes(y=cell.norm),width=box.width, outlier.shape=NA, lwd=line.width, alpha = 0.75, fill="#990000") +
          # geom_quasirandom(aes(y=cell.norm),width = 0.05, method = 'smiley',size=point.size-(0.7*point.size),alpha=0.75) + 
          # labs(title = title,subtitle ="mRNA FISH plot", y = "normalized expression",x = "genotype")
      }
    }
  }
  if(missing(plot.contrasts.2)) {
    p1 + stat_summary(aes(x=as.numeric(as.factor(genotype)) + 0.3, y=-.05),
                      fun.data = fun_length, geom = "text", size = text.size) +
      scale_x_discrete(labels = labels) + 
      theme(axis.text.x = element_text(face = "italic"))
  } else { # add secondary comparisons
    p1 + stat_summary(aes(x=genotype, y=ypos-offset), geom="text", label=plot.contrasts.2, size=text.size, colour="red") +
      stat_summary(aes(x=as.numeric(as.factor(genotype)) + 0.3, y=-.05),
                   fun.data = fun_length, geom = "text", size = text.size-1) +
      scale_x_discrete(labels = labels) + 
      theme(axis.text.x = element_text(face = "italic"))
  }
}