R/plot_tracks.R

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

#' Draw a track view of geiven genomic loci
#' @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 color named vector of colors for each sample

plot_tracks = function(bigWigs = NULL, chr, start, end, op_dir = "./", rmAfter = TRUE, bwt_path = NULL){

  check_bwtools(path = bwt_path, warn = FALSE)
  start = as.numeric(as.character(start))
  end = as.numeric(as.character(end))

  if(is.null(bigWigs)){
    stop("Please provide at least one bigWig file")
  }

  bed = data.table::data.table(chr, start, end)
  write.table(x = bed, file = "temp.bed", sep = '\t', quote = FALSE, row.names = FALSE, col.names = FALSE)

  lapply(bigWigs, function(bw){
    bn = gsub(pattern = "\\.bw$|\\.bigWig$", replacement = "", x = basename(bw), ignore.case = TRUE)
    message(paste0("Processing ", bn, ".."))

    if(!file.exists(bw)){
      stop(paste0(bw, " does not exists"))
    }
    bw = gsub(pattern = " ", replacement = "\\ ", x = bw, fixed = TRUE) #Change spaces with \ for unix style paths

    cmd = paste0("bwtool extract -locus-name jsp temp.bed ", bw, " ",  op_dir, "/temp_", bn, ".out")
    system(command = cmd, intern = TRUE)}
  )

  signals = list.files(path = op_dir, pattern = "^temp.*\\.out$", full.names = TRUE)

  signals.dat = lapply(signals, data.table::fread, sep = "\t")
  names(signals.dat) = gsub(pattern = "\\.out", replacement = "",
                            x = gsub(pattern = "temp_", replacement = "",
                                     x = basename(path = signals)))

  lapply(signals, function(x) system(paste0("rm ", x)))

  gmm = getOverlaps(chr = chr, start = start, end = end)
  gm = gmm[[1]]
  exon.tbl = gmm[[2]]

  yl = max(unlist(lapply(signals.dat, max, na.rm = TRUE)))


  par(mfrow = c(length(signals.dat)+1, 1), mgp = c(2,.33,0), las=1, tcl=-.25)

  for(i in 1:length(signals.dat)){
    x = signals.dat[[i]]
    colnames(x) = "signal"
    x[,pos := seq(start, start+nrow(x)-1)]

    mar = c(0, 3, 1, 1)
    plot(x = x$pos, y = x$signal, type = 'h', xlim = c(start, end), axes = FALSE,
         xlab = "", ylab = "", frame.plot = FALSE, ylim = c(0, yl), col = "gray70")
  }


  mar = c(0, 3, 0.5, 1)
  plot(x = 0, y = 0, xlim = c(start, end), ylim = c(0, 1),
       axes = FALSE, xlab = '', ylab = '')
  axis(side = 1, at = c(start, end), labels = c(start, end), lwd = 2,
       line = 1, cex.axis = 1.2, font = 2, lwd.ticks = 0.8)
  segments(x0 = gm[,txStart], y0 = 0.1, x1 = gm[,txEnd], y1 = 0.1, lwd = 3)
  rect(xleft = exon.tbl$start, ybottom = 0, xright = exon.tbl$end,
       ytop = 0.2, col = "gray70", border = "gray70")

}