R/plot_profile.R

In PoisonAlien/peakseason: Functions for faster ChIP/ATAC-seq analysis and visualization.

#' Draw a profile plot
#' @details This function takes output from extract_matrix and draws a profile plot
#'
#' @param mat_list Input matrix list generated by \code{\link{extract_matrix}}
#' @param summarizeBy summarize matrix by mean or median. Default mean.
#' @param condition a column name in \code{coldata} containing sample conditions. Default NULL.
#' @param collapse_reps Samples belonging to same condition to estimate average signal. Default TRUE if \code{condition} is given.
#' @param color named vector of colors for each sample
#' @param ci add SEM bars.
#' @param legend_font_size cex for font size. Default 1.
#' @param line_size line width. Default 3.
#' @param axis_lwd axis line width. Default 3.
#' @param axis_font_size font size for axis. Default 1.
#' @param sample_names names for samples. Default NULL parses from bigwig files.
#' @param auto_y_axis Default TRUE. If FALSE plots from zero to max.
#' @param y_lims Default NULL. Manually specify Y-limits.
#' @export
plot_profile = function(mat_list, summarizeBy = 'mean', color = NULL, ci = FALSE,
                        legend_font_size = 1, condition = NULL, collapse_reps = TRUE, line_size = 2, axis_lwd = 3, axis_font_size = 1,
                        sample_names = NULL, auto_y_axis = TRUE, y_lims = NULL){

  coldata = mat_list$cdata

  if(!is.null(condition)){
    if(!condition %in% colnames(coldata)){
      warning(paste0(condition, " does not exists in coldata. Here are available columns."))
      print(coldata)
      stop()
    }else{
      colnames(coldata)[which(colnames(coldata) == condition)] = "group_condition"
      coldata$group_condition = as.character(coldata$group_condition)
    }
    sample_conditons = as.character(coldata$group_condition)
  }else{
    sample_conditons = NULL
  }

  if(ci){
    ci_list = estimateCI(mats = mat_list, group = sample_conditons, collapse_reps = collapse_reps)
  }

  mat_list = summarizeMats(mats = mat_list, summarizeBy = summarizeBy, group = sample_conditons, collapse_reps = collapse_reps)

  size = as.character(mat_list$param["size"])
  mat_list = mat_list[1:(length(mat_list)-1)]

  if(is.null(condition)){
    names(mat_list) = coldata$sample_names
  }

  #return(mat_list)

  if(is.null(color)){
    color = RColorBrewer::brewer.pal(n = 8, name = "Dark2")[1:length(mat_list)]
    #names(color) = names(mat_list)
  }else{
    if(length(color) != length(mat_list)){
      warning("Insufficient number of colors. Randomly choosing few colors")
      color = c(color,
                sample(colors(), size = length(mat_list) - length(color), replace = FALSE))
    }
  }

  if(is.null(condition)){
    names(color) = names(mat_list)
  }else{
    if(collapse_reps){
      names(color) = names(mat_list)
    }else{
      cond = unique(sample_conditons)
      cond_color = color[1:length(cond)]
      names(cond_color) = cond
      print(cond_color)
      color = cond_color[sample_conditons]
    }
  }


  if(auto_y_axis){
    yl = round(c(min(unlist(lapply(mat_list, min, na.rm = TRUE))),
           max(unlist(lapply(mat_list, max, na.rm = TRUE)))), digits = 2)
  }else{
    yl = round(c(0, max(unlist(lapply(mat_list, max, na.rm = TRUE)))), digits = 2)
  }

  if(ci){
    yl[length(yl)] = yl[length(yl)] + max(unlist(lapply(ci_list, max, na.rm = TRUE)))
  }

  if(auto_y_axis){
    yl[1] = round(min(yl) - (min(yl)*0.10), digits = 2)
    yl[length(yl)] = round(max(yl) + (max(yl)*0.10), digits = 2)
  }

  if(!is.null(y_lims)){
    if(length(y_lims) < 2){
      warning("y_lim requires lower and upper limits. Auto adjusting y axis limits")
      yl[1] = round(min(yl) - (min(yl)*0.10), digits = 2)
      yl[length(yl)] = round(max(yl) + (max(yl)*0.10), digits = 2)
    }else{
      yl = sort(as.numeric(as.character(y_lims)))[1:2]
    }
  }

  xlabs = c(sapply(strsplit(x = as.character(size), split = ":", fixed = TRUE), "[[", 1), 0,
            sapply(strsplit(x = as.character(size), split = ":", fixed = TRUE), "[[", 2))
  xticks = c(0,
             as.integer(length(mat_list[[1]])/sum(as.numeric(xlabs[1]), as.numeric(xlabs[3])) * as.numeric(xlabs[1])),
             length(mat_list[[1]]))

  # mat_list = data.table::rbindlist(lapply(mat_list, function(x){
  #   x = data.table::melt(data = x)
  #   data.table::setDT(x = x)
  #   x[,xax := 1:nrow(x)]
  #   x
  # }), idcol = "sample")


  # g = ggplot(data = mat_list, aes(x = xax, y = value, col = sample))+geom_line()+
  #   scale_y_continuous(breaks = yl, labels = yl, expand = c(0, 0), limits = yl)+
  #   scale_x_continuous(breaks = xticks, labels = c(paste0("-", xlabs[1]), xlabs[2], xlabs[3]),
  #                      expand = c(0, 0), limits = xticks[c(1, 3)])+
  #   theme(panel.background  = element_blank(), axis.line.x = element_line(color="black", size = 0.2),
  #         axis.line.y = element_line(color="black", size = 0.2))
  #
  # return(g)


  par(mar = c(3, 3, 4, 2), font = 4)
  plot(mat_list[[1]], frame.plot = FALSE, axes = F,
       xlab = '', ylab = '', type = 'l',
       lwd = line_size,
       col = color[1],#color[names(mat_list)[[1]]]
       ylim = c(min(yl), max(yl)))

  if(ci){
    polygon(x = c(1:length(mat_list[[1]]), rev(1:length(mat_list[[1]]))),
            y = c(mat_list[[1]]-ci_list[[1]], rev(mat_list[[1]]+ci_list[[1]])),
            col = grDevices::adjustcolor(col = color[1], #color[names(mat_list)[[1]]],
                                         alpha.f = 0.4), border = NA)
  }
  if(length(mat_list) > 1){
    for(i in 2:length(mat_list)){

      points(mat_list[[i]], type = 'l', lwd = line_size, col = color[i])#color[names(mat_list)[[i]]])
      if(ci){
        polygon(x = c(1:length(mat_list[[i]]), rev(1:length(mat_list[[i]]))),
                y = c(mat_list[[i]]-ci_list[[i]], rev(mat_list[[i]]+ci_list[[i]])),
                col = grDevices::adjustcolor(col = color[i], #color[names(mat_list)[[i]]],
                                             alpha.f = 0.4), border = NA)
      }
    }
  }

  xlabs = c(sapply(strsplit(x = as.character(size), split = ":", fixed = TRUE), "[[", 1), 0,
            sapply(strsplit(x = as.character(size), split = ":", fixed = TRUE), "[[", 2))
  xticks = c(0,
             as.integer(length(mat_list[[1]])/sum(as.numeric(xlabs[1]), as.numeric(xlabs[3])) * as.numeric(xlabs[1])),
             length(mat_list[[1]]))

  axis(side = 1, at = xticks,
       labels = c(paste0("-", xlabs[1]), xlabs[2], xlabs[3]), lty = 1, lwd = axis_lwd,
       font = 2, cex.axis = axis_font_size, line = 0.5)
  axis(side = 2, at = yl, labels = yl, lty = 1, lwd = axis_lwd, las = 2, font = 2, cex = axis_font_size)

  if(is.null(condition)){
    add_legend(x = "topright", bty = "n", legend = names(color),
               col = color, lty = 1, lwd = 3, cex = legend_font_size, ncol = 2)
  }else{
    if(collapse_reps){
      add_legend(x = "topright", bty = "n", legend = names(color),
                 col = color, lty = 1, lwd = 3, cex = legend_font_size, ncol = 2)
    }else{
      add_legend(x = "topright", bty = "n", legend = names(cond_color),
                 col = cond_color, lty = 1, lwd = 3, cex = legend_font_size, ncol = 2)
    }
  }
}


# xpos = c(0, as.integer(mat_list[,.N,Sample][1,N]/2),  as.integer(mat_list[,.N,Sample][1,N]))


# yticks = round(seq(0, max(mat_list[, Signal], na.rm = TRUE), length.out = 4), digits = 2)

# p = ggplot(data = mat_list, aes(x = row_id, y = Signal, color = Sample))+
#   geom_line(size = 1.1)+cowplot::theme_cowplot(line_size = 1.2, font_size = 16)+
#   theme(panel.background = element_blank(), axis.text.y = element_text(size = 12), axis.text.x = element_text(size = 12))+
#   scale_x_continuous(breaks = xpos, limits = c(min(xpos), max(xpos)), labels = xlabs)+
#   scale_y_continuous(breaks = yticks, limits = c(0, max(yticks)))+
#   theme(legend.position = 'bottom', legend.title = element_blank(), legend.text = element_text(size = 10))+
#   ylab("ChIP Signal")+xlab("Position")+
#   scale_color_manual(values = color)

#  print(p)

# mat_list = data.table::rbindlist(l = lapply(mat_list, function(x) data.table::as.data.table(as.data.frame(x), keep.rownames = TRUE)),
#                                  idcol = "Sample")
# colnames(mat_list)[2:3] = c("row_id", "Signal")
# mat_list$row_id = as.numeric(as.character(gsub(pattern = "^V", replacement = "", x = mat_list$row_id)))