#' 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")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.