R/plot.deeptools.R

In robinweide/RHWlib: Simple functions to make life... simpler.

#' Plot deeptools plots.
#'
#' deeptools --outFileNameMatrix in, plot out
#'
#' @param dt A datatable from fread(mat)
#' @param zlim zlims c(min,max) for heatmap
#' @param profileZlim zlims c(min,max) for profile
#' @param subsample Sample matrix to N rows to test function
#' @param order The order to plot
#' @param ignoreSamples Remove which samples from output?
#' @param profileFunct Funciton to make profile-plots (mean, median, sum, absSum)
#' @param profileCols Vector of colors for profile-plots
#' @param sortWidth Percentage of columns to sort on (aligned on center)
#' @param sampleNames Sample names. Duh.
#' @param increasing Sort increasing?
#' @param focus Sort on which sample?
#' @param verbose Chatty?
#' @param getTopX Prints the top-X values (no plot!)
#' @param getBottomX Prints the bottom-X values  (no plot!)
#' @param heatmap Do you also want a heatmap? Tak
#' @return A plot and a df (invisible)
plot.deeptools <- function(dt, zlim, profileZlim, sampleNames, name = "", heatmap = T, smooth = F, smoothPar = 0.1, order = NULL,enrichment = F, ignoreSamples = NULL, subsample = NULL, profileCols = NULL, getTopX = NULL, getBottomX = NULL,profileFunct = "mean",sortWidth = 10, focus = 1, increasing = T, verbose = T, brewerPal = "Greys", brewerDir = 1){
  require(data.table)
  require(reshape2)
  require(ggplot2)
  require(RHWlib)
  require(dplyr)
  require(cowplot)
  # require(viridis)
  nSamples = length(sampleNames)

  dt <- as.matrix(dt)
  if(!is.null(subsample)){
    dt <- dt[sample(nrow(dt),size=subsample,replace=F),]
  }

  # which rows have no variance?


  # split mat in samples
  if(verbose){message("Splitting...")}
  perSample <- dim(dt)[2]/nSamples

  sampleList <- list()
  for(i in 1:nSamples){
    if(i == 1){
      idxes <- 1:perSample
      littleMat <- dt[,idxes]
      sampleList[[i]] <- littleMat
    } else {
      idxes <- (perSample*(i-1)+1):(perSample*(i))
      littleMat <- dt[,idxes]
      sampleList[[i]] <- littleMat
    }
  }
  names(sampleList) <- sampleNames

  # get orders beased on focus
  if(verbose){message("Sorting...")}
  foundOrder = NULL

  if(focus > nSamples){
    stop("Focus is outside of the number of samples!")
  } else {
    focusYoungGrasshopper <- sampleList[[focus]]

    # which of columns to sort on?
    nSortCol <-  (perSample/100)*sortWidth
    midpoint <- ncol(focusYoungGrasshopper)/2
    sortCols <- (midpoint-(nSortCol/2)): (midpoint+(nSortCol/2))

    foundOrder <- order(rowMeans(focusYoungGrasshopper[, sortCols]),decreasing = increasing) # in/de is flipped for ggplot
  }

  topX = NULL
  if(!is.null(getTopX)){
    topX = foundOrder[1:getTopX]
  }

  bottomX = NULL
  if(!is.null(getBottomX)){
    bottomX = foundOrder[getBottomX:length(foundOrder)]
  }

  # order individual samples
  for(i in 1:length(sampleList)){
    sampleList[[i]] <- sampleList[[i]][foundOrder,]
  }

  # melt and combining
  # Rbinding in a loop is not
  if(verbose){message("Melting...")}
  dfList <- NULL
  for(i in 1:length(sampleList)){
    if(!is.null(ignoreSamples)){
      if(i %in% ignoreSamples){next()}
    }
    melted <- reshape2::melt(t(sampleList[[i]]))
    melted$Var1 <- as.numeric(as.factor(melted$Var1))
    melted$Var2 <- as.numeric(as.factor(melted$Var2))
    melted$sample <- names(sampleList)[i]
    dfList[[i]] <- melted
  }
  df <- rbindlist(dfList)
  df$sample <- factor(df$sample, levels = sampleNames)

  # plot profiles
  groupedDF <- group_by(df, sample, Var1)

  profDF <- summarize(groupedDF, val = profileFunct(value, na.rm = T))


  if(enrichment == T){
    ZS <-  mutate(profDF, s = scale(val))
    ZS$val <- ZS$s
    profDF <- ZS[,1:3]
  }


  if(is.null(profileCols)){
    profileCols <- rep('black', nSamples)
  }

  if(!is.null(order)){
    profDF$sample <- factor(profDF$sample, levels = levels(profDF$sample)[order])
    df$sample <- factor(df$sample, levels = levels(df$sample)[order])
  }
  profDF$globalName <- name
  PRFLS <- NULL
  if(smooth == F){ # no smoothing
    PRFLS <- ggplot(profDF, aes(Var1, val, col = sample)) +
      geom_step() + facet_grid(globalName ~ sample ) + RHWtheme() +
      scale_x_continuous(expand=c(0,0), breaks = c( (perSample/2) + 0.5 ), labels = "^") +
      scale_y_continuous(expand=c(0,0), limits = profileZlim ) +
      labs( x = '', y = '') + guides(col = F) +
      scale_colour_manual(values = profileCols)
  } else {
    PRFLS <- ggplot(profDF, aes(Var1, val, col = sample)) +
      geom_smooth(method = 'loess', span = smoothPar) + facet_grid(globalName ~ sample ) + RHWtheme() +
      scale_x_continuous(expand=c(0,0), breaks = c( (perSample/2) + 0.5 ), labels = "^") +
      scale_y_continuous(expand=c(0,0), limits = profileZlim ) +
      labs( x = '', y = '') + guides(col = F) +
      scale_colour_manual(values = profileCols)
  }


  # set zlims
  if(verbose){message("Setting zlims...")}
  df[,3][df[,3] < zlim[1]] <- zlim[1]
  df[,3][df[,3] > zlim[2]] <- zlim[2]
  df$globalName <- name

  # plot heatmaps
  HTMPS <- ggplot(df, aes(Var1, Var2, fill = value)) +
    geom_raster() + facet_grid(globalName ~ sample) + RHWtheme() +
    scale_x_continuous(expand=c(0,0), breaks = c( (perSample/2) + 0.5 ), labels = "^") +
    scale_y_continuous(expand=c(0,0), breaks = NULL) +
    labs( x = '', y = '') + guides(fill = F) + scale_fill_distiller(type = 'seq', palette = brewerPal, direction = brewerDir )+
    theme(strip.text = element_blank() ,
          strip.background = element_blank(),
          plot.margin = unit( c(0,0,0,0) , units = "lines" ))

  if(heatmap){
    return(list(plot_grid(PRFLS, HTMPS, align = "v", ncol = 1, rel_heights = c(1,3)), foundOrder, profDF))
  } else {
    return(list(PRFLS, foundOrder, profDF))
  }

}