SOMNiBUS: Smooth modeling of bisulfite sequencing

Documented in binomRegMethPredPlot

#' @title Plot the predicted methylation levels
#'
#' @description This function accepts the \code{data.frame} used as an input 
#' for the function \code{binomRegMethModelPred} with additional
#' columns containing the predictions generated by the function 
#' \code{binomRegMethModelPred} and columns containing the name of each
#' experimental group and returns a plot representing the predicted methylation
#' levels according to each experimental group.
#' @param pred \code{data.frame} used as an input 
#' for the function \code{binomRegMethModelPred} (with columns 'Position', 
#' covariate names matching the original output from the 
#' function \code{binomRegMethModel}) with additional
#' columns containing the predictions generated by the function 
#' \code{binomRegMethModelPred} and columns containing the name of each
#' experimental group. Rows without a valid group name 
#' (empty character `""` or `NA`) are ignored
#' @param pred.type type of prediction returned by the function 
#' \code{binomRegMethModelPred}: \code{proportion} or \code{link.scale}.
#' The default value is "proportion".
#' @param pred.col character defines the name of the column containing the
#' prediction values. The default value is "pred". 
#' @param group.col character defines the name of the column containing the
#' experimental groups. If the group.col is set to NULL, the resulting plot 
#' will be a simple scatter plot representing all predicted values disregarding
#' any experimental design. The default value is NULL.
#' @param title the text for the title
#' @param style named list containing the wanted style 
#' (color and line type) for each experimental groups. The first
#' level list is named according each experimental group and for each 
#' experimental group there is a list containing the \code{color} and the 
#' \code{type} of the line. The line types should be among the following types:
#' \itemize{
#' \item \code{twodash},
#' \item \code{solid},
#' \item \code{longdash},
#' \item \code{dotted},
#' \item \code{dotdash},
#' \item \code{dashed},
#' \item \code{blank}.
#' }
#' The function accepts color name and its hexadecimal code. 
#' The default value is NULL meaning that the colors will be chosen randomly and
#' the line style will be set to `solid`.
#' @param save file name to create on disk. When the value is set to NULL, 
#' the plot is not saved. The default value is NULL.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information. The default value is TRUE.
#' @return This function prints out a plot of the predicted methylation levels 
#' according to preset experimental groups.
#' @author  Audrey Lemaçon
#' @importFrom data.table as.data.table
#' @importFrom tidyr pivot_longer %>%
#' @importFrom ggplot2 ggplot geom_line theme_bw aes_string theme
#' @importFrom ggplot2 geom_point geom_rug xlab ggsave ylab guides
#' @importFrom ggplot2 guide_legend scale_linetype_manual scale_color_manual
#' @examples
#' #------------------------------------------------------------#
#' data(RAdat)
#' RAdat.f <- na.omit(RAdat[RAdat$Total_Counts != 0, ])
#' BEM.obj <- binomRegMethModel(
#'   data=RAdat.f, n.k=rep(5, 3), p0=0.003307034, p1=0.9,
#'   epsilon=10^(-6), epsilon.lambda=10^(-3), maxStep=200,
#'   Quasi = FALSE, RanEff = FALSE, verbose = FALSE
#' )
#' pos <- BEM.obj$uni.pos
#' newdata <- expand.grid(pos, c(0, 1), c(0, 1))
#' colnames(newdata) <- c("Position", "T_cell", "RA")
#' my.pred <- binomRegMethModelPred(BEM.obj, newdata, type = "link.scale", 
#' verbose = FALSE)
#' newdata$group <- ""
#' newdata[(newdata$RA == 0 & newdata$T_cell == 0),]$group <- "CTRL MONO"
#' newdata[(newdata$RA == 0 & newdata$T_cell == 1),]$group <- "CTRL TCELL"
#' newdata[(newdata$RA == 1 & newdata$T_cell == 0),]$group <- "RA MONO"
#' newdata[(newdata$RA == 1 & newdata$T_cell == 1),]$group <- "RA TCELL"
#' pred <- cbind(newdata, Pred = my.pred)
#' style <- list("CTRL MONO" = list(color = "blue", type = "dashed"),
#' "CTRL TCELL" = list(color = "green", type = "dashed"),
#' "RA MONO" = list(color = "blue", type = "solid"),
#' "RA TCELL" = list(color = "green", type = "solid"))
#' g <- binomRegMethPredPlot(pred, pred.col = "Pred", group.col = "group",
#' style = style, save = NULL, verbose = FALSE)
#' 
#' @export
binomRegMethPredPlot <- function(pred, pred.type = "proportion", 
                                 pred.col = "pred", group.col = NULL, 
                                 title = "Predicted methylation levels",
                                 style = NULL, save = NULL, verbose = TRUE) {
  
  if(verbose) Message("Plot the predicted methylation levels")
  
  # check if pred.col exists in pred, raise error if not
  if(!pred.col %in% colnames(pred)){
    Error("The 'pred.col' column name is missing from the 'pred' data.frame")
  }
  
  ## test if the split list contains the 'approach' exists
  if(!pred.type %in% c("proportion","link.scale")){
    Warning(paste("Unknown prediction type. We will used",
                  "the type 'proportion'."))
    pred.type <- "proportion"
  }
  
  # if group.col == NULL, create a scatter plot with geom_point
  if(is.null(group.col)){
    g <- ggplot(data = pred, mapping = aes_string(x = "Position")) + 
      geom_point(mapping = aes_string(y = pred.col)) + 
      geom_rug(sides = "b", color = "black") + 
      theme_bw() + xlab("Genomic position") + 
      ylab(paste("Predicted methylation levels", 
                 ifelse(test = pred.type == "proportion", 
                        yes = "(proportion)", 
                        no = "(logit scale)"))) +
      ggtitle(title)
  } else {
    # check if group.col exists in pred, raise error if not
    if(!group.col %in% colnames(pred)){
      Error("The 'group.col' column name is missing from the 'pred' data.frame")
    }
    
    # suppress row missing group category
    pred <- pred[!(is.na(pred[[group.col]]) | pred[[group.col]] == ""),]
    
    # start creating the plot
    g <- ggplot(data = pred, mapping = aes_string(x = "Position")) + 
      geom_line(mapping = aes_string(y = pred.col, 
                                     group = group.col, 
                                     color = group.col,
                                     linetype = group.col
      )) + 
      geom_rug(sides = "b", color = "black") + 
      theme_bw() + xlab("Genomic position") + 
      ylab(paste("Predicted methylation levels", 
                 ifelse(test = pred.type == "proportion", 
                        yes = "(proportion)", 
                        no = "(logit scale)"))) +
      guides(color=guide_legend("Experimental groups"), 
             linetype=guide_legend("Experimental groups")) +
      ggtitle(title)
    
    # if style != NULL, check that the list names match the covs in pred,
    # raise warning if not and set style to NULL
    if(!is.null(style)){
      if(any(!unique(pred[[group.col]]) %in% names(style))){
        w_msg <- paste("Ignoring the custom style because of discrepancies",
                       "between 'style' list and the group defined",
                       "in the 'pred' data.frame.")
        Warning(w_msg)
      } else {
        style <- style[names(style) %in% unique(pred[[group.col]])]
        g <- g + scale_color_manual(breaks = names(style),
                                    values = sapply(X = style, 
                                                    FUN = function(s) s$color)) +
          scale_linetype_manual(breaks = names(style),
                                values = sapply(X = style, 
                                                FUN = function(s) s$type)) 
      }
    }
  }
  
  if(!is.null(save)){
    supported_format <- ".pdf$|png$|tiff$|bmp$|jpeg$|svg$|eps$|tex$|ps$"
    if(length(grep(supported_format, save)) == 0) {
      w_msg <- paste0("Unknown graphics format for '",save,
                      "'. The plot will be saved as a .pdf by default")
      Warning(w_msg)
      ggsave(filename = paste0(save,".pdf"), plot= g, device ="pdf")
    } else {
      ggsave(filename = save, plot = g)
    }
  }
  
  return(g)
}