R/visualizeCapRStereogene.R
In vbusa1/nearBynding: Discern RNA structure proximal to protein binding
Defines functions
visualizeCapRStereogene
Documented in
visualizeCapRStereogene
#' @title visualizeCapRStereogene
#'
#' @description Creates a visual output of all CapR RNA structure contexts
#' relative to protein binding.
#'
#' @param dir_stereogene_output Directory of stereogene output. Default working
#' directory.
#' @param CapR_prefix The prefix string common to CapR output files of
#' protein_file. Required.
#' @param protein_file A vector of at least one protein file name to be averaged
#' for visualization. File names must exclude extensions such as ".bedGraph".
#' All files in the list should be experimental or biological replicates.
#' Required.
#' @param protein_file_input A protein file name of background input to be
#' subtracted from protein_file signal. File name must exclude extension. Only
#' one input file is permitted. Optional.
#' @param x_lim A vector of two integers denoting the lower and upper x axis
#' limits. Cannot exceed wSize/2 from write_config. Default (-100, 100)
#' @param y_lim A vector of two numbers denoting the lower and upper y axis
#' limits. Optional
#' @param error A numeric value that determines the number of standard
#' deviations to show in the error bar. Default 3
#' @param nShuffle Relevant if multiple protein files are input and background
#' error has been calculated. It is the number of iterations used to derive
#' background signal error. Should be same for all protein files. Default 100.
#' @param out_file Name of output file, excluding extension. ".pdf" or ".jpeg"
#' will be added as relevant to the output file name. Default "out_file"
#' @param legend Whether a legend should be included with the output graph.
#' Default TRUE
#' @param heatmap Whether the output graph should be in the form of a heatmap
#' (TRUE) or of a line graph (FALSE). Default FALSE
#'
#' @return heatmap (JPEG) or line graph (PDF) image file
#'
#' @examples
#' ## pull example files
#' get_outfiles()
#' ## heatmap
#' visualizeCapRStereogene(CapR_prefix = "chr4and5_3UTR",
#' protein_file = "chr4and5_liftOver",
#' heatmap = TRUE,
#' out_file = "all_contexts_heatmap",
#' x_lim = c(-500, 500))
#' ## line graph
#' visualizeCapRStereogene(CapR_prefix = "chr4and5_3UTR",
#' protein_file = "chr4and5_liftOver",
#' x_lim = c(-500, 500),
#' out_file = "all_contexts_line",
#' y_lim = c(-18, 22))
#'
#' @importFrom utils read.table
#' @importFrom matrixStats rowSds
#' @importFrom graphics abline arrows lines par plot
#' @importFrom grDevices dev.off jpeg pdf rgb
#' @importFrom gplots heatmap.2 redblue
#'
#' @export
visualizeCapRStereogene <- function(dir_stereogene_output = ".",
CapR_prefix,
protein_file,
protein_file_input = NULL,
x_lim = c(-100, 100),
y_lim = NULL,
error = 1,
nShuffle = 100,
out_file = "out_file",
legend = TRUE,
heatmap = FALSE) {
if (length(protein_file) < 1) {
stop("Requires at least one protein file prefix to calculate distance")
}
if (length(protein_file) > 20) {
stop("There are > 20 protein files input. This is likely in error")
}
if (length(protein_file_input) > 1) {
stop("Only input one background track per protein.")
}
dist_bulge_1 <- NULL
for (n in seq(length(protein_file))) {
assign(paste0("dist_bulge_", n), read.table(paste0(
dir_stereogene_output,
"/", CapR_prefix, "_bulge_liftOver~",
protein_file[n], ".dist"), header = TRUE))
assign(paste0("dist_multibranch_", n), read.table(paste0(
dir_stereogene_output,
"/", CapR_prefix, "_multibranch_liftOver~",
protein_file[n], ".dist"), header = TRUE))
assign(paste0("dist_hairpin_", n), read.table(paste0(
dir_stereogene_output,
"/", CapR_prefix, "_hairpin_liftOver~",
protein_file[n], ".dist"), header = TRUE))
assign(paste0("dist_stem_", n), read.table(paste0(
dir_stereogene_output,
"/", CapR_prefix, "_stem_liftOver~",
protein_file[n], ".dist"), header = TRUE))
assign(paste0("dist_internal_", n), read.table(paste0(
dir_stereogene_output,
"/", CapR_prefix, "_internal_liftOver~",
protein_file[n], ".dist"), header = TRUE))
assign(paste0("dist_exterior_", n), read.table(paste0(
dir_stereogene_output,
"/", CapR_prefix, "_exterior_liftOver~",
protein_file[n], ".dist"), header = TRUE))
}
if (!is.null(protein_file_input)) {
dist_bulge_input <- read.table(paste0(dir_stereogene_output,
"/", CapR_prefix, "_bulge_liftOver~",
protein_file_input, ".dist"), header = TRUE)
dist_multibranch_input <- read.table(paste0(dir_stereogene_output,
"/", CapR_prefix, "_multibranch_liftOver~",
protein_file_input, ".dist"), header = TRUE)
dist_stem_input <- read.table(paste0(dir_stereogene_output,
"/", CapR_prefix, "_stem_liftOver~",
protein_file_input, ".dist"), header = TRUE)
dist_hairpin_input <- read.table(paste0(dir_stereogene_output,
"/", CapR_prefix, "_hairpin_liftOver~",
protein_file_input, ".dist"), header = TRUE)
dist_internal_input <- read.table(paste0(dir_stereogene_output,
"/", CapR_prefix, "_internal_liftOver~",
protein_file_input, ".dist"), header = TRUE)
dist_exterior_input <- read.table(paste0(dir_stereogene_output,
"/", CapR_prefix, "_exterior_liftOver~",
protein_file_input, ".dist"), header = TRUE)
}
dist_bulge <- as.data.frame(matrix(NA,
ncol = (3 * length(protein_file)) + 1,
nrow = nrow(dist_bulge_1)))
dist_stem <- as.data.frame(matrix(NA,
ncol = (3 * length(protein_file)) + 1,
nrow = nrow(dist_bulge_1)))
dist_multibranch <- as.data.frame(matrix(NA,
ncol = (3 * length(protein_file)) +1,
nrow = nrow(dist_bulge_1)))
dist_hairpin <- as.data.frame(matrix(NA,
ncol = (3 * length(protein_file)) + 1,
nrow = nrow(dist_bulge_1)))
dist_exterior <- as.data.frame(matrix(NA,
ncol = (3 * length(protein_file)) + 1,
nrow = nrow(dist_bulge_1)))
dist_internal <- as.data.frame(matrix(NA,
ncol = (3 * length(protein_file)) + 1,
nrow = nrow(dist_bulge_1)))
colnames(dist_bulge) <- c("x", paste0("Fg", seq(length(protein_file))),
paste0("Bkg", seq(length(protein_file))),
paste0("Bkg_se", seq(length(protein_file))))
colnames(dist_stem) <- c("x", paste0("Fg", seq(length(protein_file))),
paste0("Bkg", seq(length(protein_file))),
paste0("Bkg_se", seq(length(protein_file))))
colnames(dist_internal) <- c("x", paste0("Fg", seq(length(protein_file))),
paste0("Bkg", seq(length(protein_file))),
paste0("Bkg_se", seq(length(protein_file))))
colnames(dist_multibranch) <- c("x", paste0("Fg",seq(length(protein_file))),
paste0("Bkg", seq(length(protein_file))),
paste0("Bkg_se", seq(length(protein_file))))
colnames(dist_exterior) <- c("x", paste0("Fg", seq(length(protein_file))),
paste0("Bkg", seq(length(protein_file))),
paste0("Bkg_se", seq(length(protein_file))))
colnames(dist_hairpin) <- c("x", paste0("Fg", seq(length(protein_file))),
paste0("Bkg", seq(length(protein_file))),
paste0("Bkg_se", seq(length(protein_file))))
dist_hairpin$x <-dist_exterior$x <-dist_multibranch$x <-dist_internal$x <-
dist_stem$x <-dist_bulge$x <-dist_bulge_1$x
for (n in seq(length(protein_file))) {
dist_bulge[, 1 + n] <-
eval(parse(text = paste0("dist_bulge_", n)))$Fg
dist_stem[, 1 + n] <-
eval(parse(text = paste0("dist_stem_", n)))$Fg
dist_hairpin[, 1 + n] <-
eval(parse(text = paste0("dist_hairpin_", n)))$Fg
dist_multibranch[, 1 + n] <-
eval(parse(text = paste0("dist_multibranch_", n)))$Fg
dist_exterior[, 1 + n] <-
eval(parse(text = paste0("dist_exterior_", n)))$Fg
dist_internal[, 1 + n] <-
eval(parse(text = paste0("dist_internal_", n)))$Fg
dist_bulge[, 1 + n + length(protein_file)] <-
eval(parse(text = paste0("dist_bulge_", n)))$Bkg
dist_stem[, 1 + n + length(protein_file)] <-
eval(parse(text = paste0("dist_stem_", n)))$Bkg
dist_hairpin[, 1 + n + length(protein_file)] <-
eval(parse(text = paste0("dist_hairpin_", n)))$Bkg
dist_multibranch[, 1 + n + length(protein_file)] <-
eval(parse(text = paste0("dist_multibranch_", n)))$Bkg
dist_exterior[, 1 + n + length(protein_file)] <-
eval(parse(text = paste0("dist_exterior_", n)))$Bkg
dist_internal[, 1 + n + length(protein_file)] <-
eval(parse(text = paste0("dist_internal_", n)))$Bkg
dist_bulge[, 1 + n + 2*length(protein_file)] <-
eval(parse(text = paste0("dist_bulge_", n)))$Bkg_se
dist_stem[, 1 + n + 2*length(protein_file)] <-
eval(parse(text = paste0("dist_stem_", n)))$Bkg_se
dist_hairpin[, 1 + n + 2*length(protein_file)] <-
eval(parse(text = paste0("dist_hairpin_", n)))$Bkg_se
dist_multibranch[, 1 + n + 2*length(protein_file)] <-
eval(parse(text = paste0("dist_multibranch_", n)))$Bkg_se
dist_exterior[, 1 + n + 2*length(protein_file)] <-
eval(parse(text = paste0("dist_exterior_", n)))$Bkg_se
dist_internal[, 1 + n + 2*length(protein_file)] <-
eval(parse(text = paste0("dist_internal_", n)))$Bkg_se
}
if (!is.null(protein_file_input)) {
dist_bulge[, 2:(length(protein_file) + 1)] <-
dist_bulge[, 2:(length(protein_file) + 1)] -
dist_bulge_input$Fg
dist_bulge[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] <-
dist_bulge[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] -
dist_bulge_input$Bkg
dist_stem[,2:(length(protein_file) + 1)] <-
dist_stem[, 2:(length(protein_file) + 1)] -
dist_stem_input$Fg
dist_stem[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] <-
dist_stem[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] -
dist_stem_input$Bkg
dist_hairpin[, 2:(length(protein_file) + 1)] <-
dist_hairpin[, 2:(length(protein_file) + 1)] -
dist_hairpin_input$Fg
dist_hairpin[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] <-
dist_hairpin[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] -
dist_hairpin_input$Bkg
dist_multibranch[, 2:(length(protein_file) + 1)] <-
dist_multibranch[,2:(length(protein_file) + 1)] -
dist_multibranch_input$Fg
dist_multibranch[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] <-
dist_multibranch[,
(length(protein_file)+2):((length(protein_file) * 2) + 1)] -
dist_multibranch_input$Bkg
dist_internal[, 2:(length(protein_file) + 1)] <-
dist_internal[, 2:(length(protein_file) + 1)] -
dist_internal_input$Fg
dist_internal[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] <-
dist_internal[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] -
dist_internal_input$Bkg
dist_exterior[, 2:(length(protein_file) + 1)] <-
dist_exterior[, 2:(length(protein_file) + 1)] -
dist_exterior_input$Fg
dist_exterior[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] <-
dist_exterior[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] -
dist_exterior_input$Bkg
}
if (length(protein_file) == 1) {
dist_bulge$Fg <- dist_bulge$Fg1
dist_bulge$Fg_se <- 0
dist_bulge$Bkg <- dist_bulge$Bkg1
dist_bulge$Bkg_se <- dist_bulge$Bkg_se1
dist_hairpin$Fg <- dist_hairpin$Fg1
dist_hairpin$Fg_se <- 0
dist_hairpin$Bkg <- dist_hairpin$Bkg1
dist_hairpin$Bkg_se <- dist_hairpin$Bkg_se1
dist_internal$Fg <- dist_internal$Fg1
dist_internal$Fg_se <- 0
dist_internal$Bkg <- dist_internal$Bkg1
dist_internal$Bkg_se <- dist_internal$Bkg_se1
dist_exterior$Fg <- dist_exterior$Fg1
dist_exterior$Fg_se <- 0
dist_exterior$Bkg <- dist_exterior$Bkg1
dist_exterior$Bkg_se <- dist_exterior$Bkg_se1
dist_stem$Fg <- dist_stem$Fg1
dist_stem$Fg_se <- 0
dist_stem$Bkg <- dist_stem$Bkg1
dist_stem$Bkg_se <- dist_stem$Bkg_se1
dist_multibranch$Fg <- dist_multibranch$Fg1
dist_multibranch$Fg_se <- 0
dist_multibranch$Bkg <- dist_multibranch$Bkg1
dist_multibranch$Bkg_se <- dist_multibranch$Bkg_se1
} else {
dist_bulge$Fg <- rowMeans(dist_bulge[, 2:(length(protein_file) + 1)])
dist_bulge$Fg_se <- rowSds(as.matrix(dist_bulge[,
2:(length(protein_file) + 1)]))/sqrt(length(protein_file))
dist_bulge$Bkg <- rowMeans(dist_bulge[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)])
dist_bulge$Bkg_se <- apply(as.matrix(dist_bulge[,
((length(protein_file) * 2) + 2):((length(protein_file) * 3) + 1)]),
1, function(x){
x<-x*sqrt(nShuffle)
return(sqrt(sum(x*x)/nShuffle))}) +
rowSds(as.matrix(dist_bulge[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)]))/
sqrt(length(protein_file))
dist_hairpin$Fg <- rowMeans(dist_hairpin[,
2:(length(protein_file) + 1)])
dist_hairpin$Fg_se <- rowSds(as.matrix(dist_hairpin[,
2:(length(protein_file) + 1)]))/sqrt(length(protein_file))
dist_hairpin$Bkg <- rowMeans(dist_hairpin[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)])
dist_hairpin$Bkg_se <- apply(as.matrix(dist_hairpin[,
((length(protein_file) * 2) + 2):((length(protein_file) * 3) + 1)]),
1, function(x){
x<-x*sqrt(nShuffle)
return(sqrt(sum(x*x)/nShuffle))}) +
rowSds(as.matrix(dist_hairpin[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)]))/
sqrt(length(protein_file))
dist_stem$Fg <- rowMeans(dist_stem[, 2:(length(protein_file) + 1)])
dist_stem$Fg_se <- rowSds(as.matrix(dist_stem[,
2:(length(protein_file) + 1)]))/sqrt(length(protein_file))
dist_stem$Bkg <- rowMeans(dist_stem[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)])
dist_stem$Bkg_se <- apply(as.matrix(dist_stem[,
((length(protein_file) * 2) + 2):((length(protein_file) * 3) + 1)]),
1, function(x){
x<-x*sqrt(nShuffle)
return(sqrt(sum(x*x)/nShuffle))}) +
rowSds(as.matrix(dist_stem[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)]))/
sqrt(length(protein_file))
dist_multibranch$Fg<-rowMeans(dist_multibranch[,
2:(length(protein_file)+ 1)])
dist_multibranch$Fg_se <- rowSds(as.matrix(dist_multibranch[,
2:(length(protein_file) + 1)]))/sqrt(length(protein_file))
dist_multibranch$Bkg <- rowMeans(dist_multibranch[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)])
dist_multibranch$Bkg_se <- apply(as.matrix(dist_multibranch[,
((length(protein_file) * 2) + 2):((length(protein_file) * 3) + 1)]),
1, function(x){
x<-x*sqrt(nShuffle)
return(sqrt(sum(x*x)/nShuffle))}) +
rowSds(as.matrix(dist_multibranch[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)]))/
sqrt(length(protein_file))
dist_internal$Fg <- rowMeans(dist_internal[,
2:(length(protein_file) + 1)])
dist_internal$Fg_se <- rowSds(as.matrix(dist_internal[,
2:(length(protein_file) + 1)]))/sqrt(length(protein_file))
dist_internal$Bkg <- rowMeans(dist_internal[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)])
dist_internal$Bkg_se <- apply(as.matrix(dist_internal[,
((length(protein_file) * 2) + 2):((length(protein_file) * 3) + 1)]),
1, function(x){
x<-x*sqrt(nShuffle)
return(sqrt(sum(x*x)/nShuffle))}) +
rowSds(as.matrix(dist_internal[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)]))/
sqrt(length(protein_file))
dist_exterior$Fg <- rowMeans(dist_exterior[,
2:(length(protein_file) + 1)])
dist_exterior$Fg_se <- rowSds(as.matrix(dist_exterior[,
2:(length(protein_file) + 1)]))/sqrt(length(protein_file))
dist_exterior$Bkg <- rowMeans(dist_exterior[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)])
dist_exterior$Bkg_se <- apply(as.matrix(dist_exterior[,
((length(protein_file) * 2) + 2):((length(protein_file) * 3) + 1)]),
1, function(x){
x<-x*sqrt(nShuffle)
return(sqrt(sum(x*x)/nShuffle))}) +
rowSds(as.matrix(dist_exterior[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)]))/
sqrt(length(protein_file))
}
Bkg<-data.frame(Bkg = rowMeans(cbind(dist_multibranch$Bkg,
dist_stem$Bkg,
dist_bulge$Bkg,
dist_hairpin$Bkg,
dist_exterior$Bkg,
dist_internal$Bkg)),
Bkg_se = pmax(dist_multibranch$Bkg_se,
dist_stem$Bkg_se,
dist_bulge$Bkg_se,
dist_hairpin$Bkg_se,
dist_exterior$Bkg_se,
dist_internal$Bkg_se))
if (heatmap == FALSE) {
# create line plot
out_file <- paste0(out_file, ".pdf")
# save plot to pdf
pdf(out_file, height = 4, width = 6)
# create the plot
old.par <- par(no.readonly = TRUE)
par(oma = c(0, 0, 0, 0), mar = c(3, 3, 3, 1), mgp = c(1.6, 0.45, 0))
if (is.null(y_lim)) {
max_y <- max(c(max(dist_exterior$Fg), max(dist_internal$Fg),
max(dist_stem$Fg), max(dist_hairpin$Fg),
max(dist_multibranch$Fg), max(dist_bulge$Fg))) + 0.1
min_y <- min(c(min(dist_exterior$Fg), min(dist_internal$Fg),
min(dist_stem$Fg), min(dist_hairpin$Fg),
min(dist_multibranch$Fg), min(dist_bulge$Fg))) - 0.1
y_lim <- c(min_y, max_y)
}
plot(dist_bulge$x, Bkg$Bkg, type = "l",
col = "black", xlim = x_lim, ylim = y_lim,
main = NULL, xlab = "Nucleotide", ylab = "Density x 100",
cex.axis = 0.8, cex.lab = 1, cex.main = 1.2, lwd = 2)
arrows(dist_bulge$x, Bkg$Bkg - (error * Bkg$Bkg_se),
dist_bulge$x, Bkg$Bkg + (error * Bkg$Bkg_se),
length = 0, angle = 90, code = 3,
col = rgb(blue = 0, red = 0, green = 0, alpha = 0.5))
abline(v = 0, col = "grey", lty = 2)
lines(dist_bulge$x, dist_bulge$Fg, col = "blue", lwd = 2)
lines(dist_multibranch$x, dist_multibranch$Fg, col = "red", lwd = 2)
lines(dist_stem$x, dist_stem$Fg, col = "green", lwd = 2)
lines(dist_hairpin$x, dist_hairpin$Fg, col = "purple", lwd = 2)
lines(dist_internal$x, dist_internal$Fg, col = "orange", lwd = 2)
lines(dist_exterior$x, dist_exterior$Fg, col = "cyan", lwd = 2)
arrows(dist_bulge$x, dist_bulge$Fg - (error * dist_bulge$Fg_se),
dist_bulge$x, dist_bulge$Fg + (error * dist_bulge$Fg_se),
length = 0, angle = 90, code = 3,
col = rgb(blue = 1, red = 0, green = 0, alpha = 0.5))
arrows(dist_multibranch$x, dist_multibranch$Fg -
(error * dist_multibranch$Fg_se), dist_multibranch$x,
dist_multibranch$Fg + (error * dist_multibranch$Fg_se),
length = 0, angle = 90, code = 3,
col = rgb(blue = 0, red = 1, green = 0, alpha = 0.5))
arrows(dist_stem$x, dist_stem$Fg - (error * dist_stem$Fg_se),
dist_stem$x, dist_stem$Fg + (error * dist_stem$Fg_se),
length = 0, angle = 90, code = 3,
col = rgb(blue = 0, red = 0, green = 1, alpha = 0.5))
arrows(dist_hairpin$x, dist_hairpin$Fg - (error * dist_hairpin$Fg_se),
dist_hairpin$x, dist_hairpin$Fg + (error * dist_hairpin$Fg_se),
length = 0, angle = 90, code = 3,
col = rgb(blue = 1, red = 1, green = 0, alpha = 0.5))
arrows(dist_internal$x, dist_internal$Fg -
(error*dist_internal$Fg_se), dist_internal$x, dist_internal$Fg +
(error * dist_internal$Fg_se), length = 0, angle = 90,
code = 3, col = rgb(blue = 0, red = 1, green = 1, alpha = 0.5))
arrows(dist_exterior$x, dist_exterior$Fg -
(error * dist_exterior$Fg_se), dist_exterior$x,
dist_exterior$Fg +
(error * dist_exterior$Fg_se), length = 0, angle = 90,
code = 3, col = rgb(blue = 1, red = 0, green = 1, alpha = 0.5))
if (legend == TRUE) {
legend("bottomright",
legend = c("Background", "Bulge", "Multibranch", "Stem",
"Hairpin", "Internal", "Exterior"),
col = c("black", "blue", "red", "green", "purple",
"orange", "cyan"), lty = 1, cex = 0.8)
}
dev.off()
} else {
out_file <- paste0(out_file, ".jpeg")
coords_1 <- which(dist_stem$x == x_lim[1])
coords_2 <- which(dist_stem$x == x_lim[2])
heatmap_plot <- data.frame(stem = dist_stem$Fg[coords_1:coords_2],
hairpin = dist_hairpin$Fg[coords_1:coords_2],
bulge = dist_bulge$Fg[coords_1:coords_2],
internal = dist_internal$Fg[coords_1:coords_2],
multibranch = dist_multibranch$Fg[coords_1:coords_2],
exterior = dist_exterior$Fg[coords_1:coords_2])
heatmap_plot <- as.matrix(t(heatmap_plot))
zero <- which(dist_stem$x == 0) - coords_1
if(zero < 0){
colsep <- 0
} else{
colsep <- c(zero - 1, zero)
}
if (is.null(y_lim)) {
max_y <- max(c(max(dist_exterior$Fg), max(dist_internal$Fg),
max(dist_stem$Fg), max(dist_hairpin$Fg),
max(dist_multibranch$Fg), max(dist_bulge$Fg))) + 0.1
min_y <- min(c(min(dist_exterior$Fg), min(dist_internal$Fg),
min(dist_stem$Fg), min(dist_hairpin$Fg),
min(dist_multibranch$Fg), min(dist_bulge$Fg))) - 0.1
y_lim <- c(-max(abs(max_y), abs(min_y)),
max(abs(max_y), abs(min_y)))
}
key <- TRUE
if (legend == FALSE) {key <- FALSE}
# make plot
jpeg(out_file, height = 4, width = 8, units = "in", res = 300)
heatmap.2(heatmap_plot, Rowv = NA, Colv = NA,
col = redblue(256), dendrogram = "none",
labCol = FALSE, trace = "none", density.info = "none",
symkey = FALSE, key = key, colsep = colsep,
sepcolor = "grey", margins = c(1, 8),
breaks = seq(y_lim[1], y_lim[2], length.out = 257))
dev.off()
}
}
vbusa1/nearBynding documentation built on Aug. 4, 2021, 4:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions visualizeCapRStereogene
Documented in visualizeCapRStereogene
#' @title visualizeCapRStereogene
#'
#' @description Creates a visual output of all CapR RNA structure contexts
#' relative to protein binding.
#'
#' @param dir_stereogene_output Directory of stereogene output. Default working
#' directory.
#' @param CapR_prefix The prefix string common to CapR output files of
#' protein_file. Required.
#' @param protein_file A vector of at least one protein file name to be averaged
#' for visualization. File names must exclude extensions such as ".bedGraph".
#' All files in the list should be experimental or biological replicates.
#' Required.
#' @param protein_file_input A protein file name of background input to be
#' subtracted from protein_file signal. File name must exclude extension. Only
#' one input file is permitted. Optional.
#' @param x_lim A vector of two integers denoting the lower and upper x axis
#' limits. Cannot exceed wSize/2 from write_config. Default (-100, 100)
#' @param y_lim A vector of two numbers denoting the lower and upper y axis
#' limits. Optional
#' @param error A numeric value that determines the number of standard
#' deviations to show in the error bar. Default 3
#' @param nShuffle Relevant if multiple protein files are input and background
#' error has been calculated. It is the number of iterations used to derive
#' background signal error. Should be same for all protein files. Default 100.
#' @param out_file Name of output file, excluding extension. ".pdf" or ".jpeg"
#' will be added as relevant to the output file name. Default "out_file"
#' @param legend Whether a legend should be included with the output graph.
#' Default TRUE
#' @param heatmap Whether the output graph should be in the form of a heatmap
#' (TRUE) or of a line graph (FALSE). Default FALSE
#'
#' @return heatmap (JPEG) or line graph (PDF) image file
#'
#' @examples
#' ## pull example files
#' get_outfiles()
#' ## heatmap
#' visualizeCapRStereogene(CapR_prefix = "chr4and5_3UTR",
#' protein_file = "chr4and5_liftOver",
#' heatmap = TRUE,
#' out_file = "all_contexts_heatmap",
#' x_lim = c(-500, 500))
#' ## line graph
#' visualizeCapRStereogene(CapR_prefix = "chr4and5_3UTR",
#' protein_file = "chr4and5_liftOver",
#' x_lim = c(-500, 500),
#' out_file = "all_contexts_line",
#' y_lim = c(-18, 22))
#'
#' @importFrom utils read.table
#' @importFrom matrixStats rowSds
#' @importFrom graphics abline arrows lines par plot
#' @importFrom grDevices dev.off jpeg pdf rgb
#' @importFrom gplots heatmap.2 redblue
#'
#' @export
visualizeCapRStereogene <- function(dir_stereogene_output = ".",
CapR_prefix,
protein_file,
protein_file_input = NULL,
x_lim = c(-100, 100),
y_lim = NULL,
error = 1,
nShuffle = 100,
out_file = "out_file",
legend = TRUE,
heatmap = FALSE) {
if (length(protein_file) < 1) {
stop("Requires at least one protein file prefix to calculate distance")
}
if (length(protein_file) > 20) {
stop("There are > 20 protein files input. This is likely in error")
}
if (length(protein_file_input) > 1) {
stop("Only input one background track per protein.")
}
dist_bulge_1 <- NULL
for (n in seq(length(protein_file))) {
assign(paste0("dist_bulge_", n), read.table(paste0(
dir_stereogene_output,
"/", CapR_prefix, "_bulge_liftOver~",
protein_file[n], ".dist"), header = TRUE))
assign(paste0("dist_multibranch_", n), read.table(paste0(
dir_stereogene_output,
"/", CapR_prefix, "_multibranch_liftOver~",
protein_file[n], ".dist"), header = TRUE))
assign(paste0("dist_hairpin_", n), read.table(paste0(
dir_stereogene_output,
"/", CapR_prefix, "_hairpin_liftOver~",
protein_file[n], ".dist"), header = TRUE))
assign(paste0("dist_stem_", n), read.table(paste0(
dir_stereogene_output,
"/", CapR_prefix, "_stem_liftOver~",
protein_file[n], ".dist"), header = TRUE))
assign(paste0("dist_internal_", n), read.table(paste0(
dir_stereogene_output,
"/", CapR_prefix, "_internal_liftOver~",
protein_file[n], ".dist"), header = TRUE))
assign(paste0("dist_exterior_", n), read.table(paste0(
dir_stereogene_output,
"/", CapR_prefix, "_exterior_liftOver~",
protein_file[n], ".dist"), header = TRUE))
}
if (!is.null(protein_file_input)) {
dist_bulge_input <- read.table(paste0(dir_stereogene_output,
"/", CapR_prefix, "_bulge_liftOver~",
protein_file_input, ".dist"), header = TRUE)
dist_multibranch_input <- read.table(paste0(dir_stereogene_output,
"/", CapR_prefix, "_multibranch_liftOver~",
protein_file_input, ".dist"), header = TRUE)
dist_stem_input <- read.table(paste0(dir_stereogene_output,
"/", CapR_prefix, "_stem_liftOver~",
protein_file_input, ".dist"), header = TRUE)
dist_hairpin_input <- read.table(paste0(dir_stereogene_output,
"/", CapR_prefix, "_hairpin_liftOver~",
protein_file_input, ".dist"), header = TRUE)
dist_internal_input <- read.table(paste0(dir_stereogene_output,
"/", CapR_prefix, "_internal_liftOver~",
protein_file_input, ".dist"), header = TRUE)
dist_exterior_input <- read.table(paste0(dir_stereogene_output,
"/", CapR_prefix, "_exterior_liftOver~",
protein_file_input, ".dist"), header = TRUE)
}
dist_bulge <- as.data.frame(matrix(NA,
ncol = (3 * length(protein_file)) + 1,
nrow = nrow(dist_bulge_1)))
dist_stem <- as.data.frame(matrix(NA,
ncol = (3 * length(protein_file)) + 1,
nrow = nrow(dist_bulge_1)))
dist_multibranch <- as.data.frame(matrix(NA,
ncol = (3 * length(protein_file)) +1,
nrow = nrow(dist_bulge_1)))
dist_hairpin <- as.data.frame(matrix(NA,
ncol = (3 * length(protein_file)) + 1,
nrow = nrow(dist_bulge_1)))
dist_exterior <- as.data.frame(matrix(NA,
ncol = (3 * length(protein_file)) + 1,
nrow = nrow(dist_bulge_1)))
dist_internal <- as.data.frame(matrix(NA,
ncol = (3 * length(protein_file)) + 1,
nrow = nrow(dist_bulge_1)))
colnames(dist_bulge) <- c("x", paste0("Fg", seq(length(protein_file))),
paste0("Bkg", seq(length(protein_file))),
paste0("Bkg_se", seq(length(protein_file))))
colnames(dist_stem) <- c("x", paste0("Fg", seq(length(protein_file))),
paste0("Bkg", seq(length(protein_file))),
paste0("Bkg_se", seq(length(protein_file))))
colnames(dist_internal) <- c("x", paste0("Fg", seq(length(protein_file))),
paste0("Bkg", seq(length(protein_file))),
paste0("Bkg_se", seq(length(protein_file))))
colnames(dist_multibranch) <- c("x", paste0("Fg",seq(length(protein_file))),
paste0("Bkg", seq(length(protein_file))),
paste0("Bkg_se", seq(length(protein_file))))
colnames(dist_exterior) <- c("x", paste0("Fg", seq(length(protein_file))),
paste0("Bkg", seq(length(protein_file))),
paste0("Bkg_se", seq(length(protein_file))))
colnames(dist_hairpin) <- c("x", paste0("Fg", seq(length(protein_file))),
paste0("Bkg", seq(length(protein_file))),
paste0("Bkg_se", seq(length(protein_file))))
dist_hairpin$x <-dist_exterior$x <-dist_multibranch$x <-dist_internal$x <-
dist_stem$x <-dist_bulge$x <-dist_bulge_1$x
for (n in seq(length(protein_file))) {
dist_bulge[, 1 + n] <-
eval(parse(text = paste0("dist_bulge_", n)))$Fg
dist_stem[, 1 + n] <-
eval(parse(text = paste0("dist_stem_", n)))$Fg
dist_hairpin[, 1 + n] <-
eval(parse(text = paste0("dist_hairpin_", n)))$Fg
dist_multibranch[, 1 + n] <-
eval(parse(text = paste0("dist_multibranch_", n)))$Fg
dist_exterior[, 1 + n] <-
eval(parse(text = paste0("dist_exterior_", n)))$Fg
dist_internal[, 1 + n] <-
eval(parse(text = paste0("dist_internal_", n)))$Fg
dist_bulge[, 1 + n + length(protein_file)] <-
eval(parse(text = paste0("dist_bulge_", n)))$Bkg
dist_stem[, 1 + n + length(protein_file)] <-
eval(parse(text = paste0("dist_stem_", n)))$Bkg
dist_hairpin[, 1 + n + length(protein_file)] <-
eval(parse(text = paste0("dist_hairpin_", n)))$Bkg
dist_multibranch[, 1 + n + length(protein_file)] <-
eval(parse(text = paste0("dist_multibranch_", n)))$Bkg
dist_exterior[, 1 + n + length(protein_file)] <-
eval(parse(text = paste0("dist_exterior_", n)))$Bkg
dist_internal[, 1 + n + length(protein_file)] <-
eval(parse(text = paste0("dist_internal_", n)))$Bkg
dist_bulge[, 1 + n + 2*length(protein_file)] <-
eval(parse(text = paste0("dist_bulge_", n)))$Bkg_se
dist_stem[, 1 + n + 2*length(protein_file)] <-
eval(parse(text = paste0("dist_stem_", n)))$Bkg_se
dist_hairpin[, 1 + n + 2*length(protein_file)] <-
eval(parse(text = paste0("dist_hairpin_", n)))$Bkg_se
dist_multibranch[, 1 + n + 2*length(protein_file)] <-
eval(parse(text = paste0("dist_multibranch_", n)))$Bkg_se
dist_exterior[, 1 + n + 2*length(protein_file)] <-
eval(parse(text = paste0("dist_exterior_", n)))$Bkg_se
dist_internal[, 1 + n + 2*length(protein_file)] <-
eval(parse(text = paste0("dist_internal_", n)))$Bkg_se
}
if (!is.null(protein_file_input)) {
dist_bulge[, 2:(length(protein_file) + 1)] <-
dist_bulge[, 2:(length(protein_file) + 1)] -
dist_bulge_input$Fg
dist_bulge[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] <-
dist_bulge[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] -
dist_bulge_input$Bkg
dist_stem[,2:(length(protein_file) + 1)] <-
dist_stem[, 2:(length(protein_file) + 1)] -
dist_stem_input$Fg
dist_stem[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] <-
dist_stem[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] -
dist_stem_input$Bkg
dist_hairpin[, 2:(length(protein_file) + 1)] <-
dist_hairpin[, 2:(length(protein_file) + 1)] -
dist_hairpin_input$Fg
dist_hairpin[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] <-
dist_hairpin[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] -
dist_hairpin_input$Bkg
dist_multibranch[, 2:(length(protein_file) + 1)] <-
dist_multibranch[,2:(length(protein_file) + 1)] -
dist_multibranch_input$Fg
dist_multibranch[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] <-
dist_multibranch[,
(length(protein_file)+2):((length(protein_file) * 2) + 1)] -
dist_multibranch_input$Bkg
dist_internal[, 2:(length(protein_file) + 1)] <-
dist_internal[, 2:(length(protein_file) + 1)] -
dist_internal_input$Fg
dist_internal[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] <-
dist_internal[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] -
dist_internal_input$Bkg
dist_exterior[, 2:(length(protein_file) + 1)] <-
dist_exterior[, 2:(length(protein_file) + 1)] -
dist_exterior_input$Fg
dist_exterior[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] <-
dist_exterior[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)] -
dist_exterior_input$Bkg
}
if (length(protein_file) == 1) {
dist_bulge$Fg <- dist_bulge$Fg1
dist_bulge$Fg_se <- 0
dist_bulge$Bkg <- dist_bulge$Bkg1
dist_bulge$Bkg_se <- dist_bulge$Bkg_se1
dist_hairpin$Fg <- dist_hairpin$Fg1
dist_hairpin$Fg_se <- 0
dist_hairpin$Bkg <- dist_hairpin$Bkg1
dist_hairpin$Bkg_se <- dist_hairpin$Bkg_se1
dist_internal$Fg <- dist_internal$Fg1
dist_internal$Fg_se <- 0
dist_internal$Bkg <- dist_internal$Bkg1
dist_internal$Bkg_se <- dist_internal$Bkg_se1
dist_exterior$Fg <- dist_exterior$Fg1
dist_exterior$Fg_se <- 0
dist_exterior$Bkg <- dist_exterior$Bkg1
dist_exterior$Bkg_se <- dist_exterior$Bkg_se1
dist_stem$Fg <- dist_stem$Fg1
dist_stem$Fg_se <- 0
dist_stem$Bkg <- dist_stem$Bkg1
dist_stem$Bkg_se <- dist_stem$Bkg_se1
dist_multibranch$Fg <- dist_multibranch$Fg1
dist_multibranch$Fg_se <- 0
dist_multibranch$Bkg <- dist_multibranch$Bkg1
dist_multibranch$Bkg_se <- dist_multibranch$Bkg_se1
} else {
dist_bulge$Fg <- rowMeans(dist_bulge[, 2:(length(protein_file) + 1)])
dist_bulge$Fg_se <- rowSds(as.matrix(dist_bulge[,
2:(length(protein_file) + 1)]))/sqrt(length(protein_file))
dist_bulge$Bkg <- rowMeans(dist_bulge[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)])
dist_bulge$Bkg_se <- apply(as.matrix(dist_bulge[,
((length(protein_file) * 2) + 2):((length(protein_file) * 3) + 1)]),
1, function(x){
x<-x*sqrt(nShuffle)
return(sqrt(sum(x*x)/nShuffle))}) +
rowSds(as.matrix(dist_bulge[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)]))/
sqrt(length(protein_file))
dist_hairpin$Fg <- rowMeans(dist_hairpin[,
2:(length(protein_file) + 1)])
dist_hairpin$Fg_se <- rowSds(as.matrix(dist_hairpin[,
2:(length(protein_file) + 1)]))/sqrt(length(protein_file))
dist_hairpin$Bkg <- rowMeans(dist_hairpin[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)])
dist_hairpin$Bkg_se <- apply(as.matrix(dist_hairpin[,
((length(protein_file) * 2) + 2):((length(protein_file) * 3) + 1)]),
1, function(x){
x<-x*sqrt(nShuffle)
return(sqrt(sum(x*x)/nShuffle))}) +
rowSds(as.matrix(dist_hairpin[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)]))/
sqrt(length(protein_file))
dist_stem$Fg <- rowMeans(dist_stem[, 2:(length(protein_file) + 1)])
dist_stem$Fg_se <- rowSds(as.matrix(dist_stem[,
2:(length(protein_file) + 1)]))/sqrt(length(protein_file))
dist_stem$Bkg <- rowMeans(dist_stem[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)])
dist_stem$Bkg_se <- apply(as.matrix(dist_stem[,
((length(protein_file) * 2) + 2):((length(protein_file) * 3) + 1)]),
1, function(x){
x<-x*sqrt(nShuffle)
return(sqrt(sum(x*x)/nShuffle))}) +
rowSds(as.matrix(dist_stem[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)]))/
sqrt(length(protein_file))
dist_multibranch$Fg<-rowMeans(dist_multibranch[,
2:(length(protein_file)+ 1)])
dist_multibranch$Fg_se <- rowSds(as.matrix(dist_multibranch[,
2:(length(protein_file) + 1)]))/sqrt(length(protein_file))
dist_multibranch$Bkg <- rowMeans(dist_multibranch[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)])
dist_multibranch$Bkg_se <- apply(as.matrix(dist_multibranch[,
((length(protein_file) * 2) + 2):((length(protein_file) * 3) + 1)]),
1, function(x){
x<-x*sqrt(nShuffle)
return(sqrt(sum(x*x)/nShuffle))}) +
rowSds(as.matrix(dist_multibranch[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)]))/
sqrt(length(protein_file))
dist_internal$Fg <- rowMeans(dist_internal[,
2:(length(protein_file) + 1)])
dist_internal$Fg_se <- rowSds(as.matrix(dist_internal[,
2:(length(protein_file) + 1)]))/sqrt(length(protein_file))
dist_internal$Bkg <- rowMeans(dist_internal[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)])
dist_internal$Bkg_se <- apply(as.matrix(dist_internal[,
((length(protein_file) * 2) + 2):((length(protein_file) * 3) + 1)]),
1, function(x){
x<-x*sqrt(nShuffle)
return(sqrt(sum(x*x)/nShuffle))}) +
rowSds(as.matrix(dist_internal[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)]))/
sqrt(length(protein_file))
dist_exterior$Fg <- rowMeans(dist_exterior[,
2:(length(protein_file) + 1)])
dist_exterior$Fg_se <- rowSds(as.matrix(dist_exterior[,
2:(length(protein_file) + 1)]))/sqrt(length(protein_file))
dist_exterior$Bkg <- rowMeans(dist_exterior[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)])
dist_exterior$Bkg_se <- apply(as.matrix(dist_exterior[,
((length(protein_file) * 2) + 2):((length(protein_file) * 3) + 1)]),
1, function(x){
x<-x*sqrt(nShuffle)
return(sqrt(sum(x*x)/nShuffle))}) +
rowSds(as.matrix(dist_exterior[,
(length(protein_file) + 2):((length(protein_file) * 2) + 1)]))/
sqrt(length(protein_file))
}
Bkg<-data.frame(Bkg = rowMeans(cbind(dist_multibranch$Bkg,
dist_stem$Bkg,
dist_bulge$Bkg,
dist_hairpin$Bkg,
dist_exterior$Bkg,
dist_internal$Bkg)),
Bkg_se = pmax(dist_multibranch$Bkg_se,
dist_stem$Bkg_se,
dist_bulge$Bkg_se,
dist_hairpin$Bkg_se,
dist_exterior$Bkg_se,
dist_internal$Bkg_se))
if (heatmap == FALSE) {
# create line plot
out_file <- paste0(out_file, ".pdf")
# save plot to pdf
pdf(out_file, height = 4, width = 6)
# create the plot
old.par <- par(no.readonly = TRUE)
par(oma = c(0, 0, 0, 0), mar = c(3, 3, 3, 1), mgp = c(1.6, 0.45, 0))
if (is.null(y_lim)) {
max_y <- max(c(max(dist_exterior$Fg), max(dist_internal$Fg),
max(dist_stem$Fg), max(dist_hairpin$Fg),
max(dist_multibranch$Fg), max(dist_bulge$Fg))) + 0.1
min_y <- min(c(min(dist_exterior$Fg), min(dist_internal$Fg),
min(dist_stem$Fg), min(dist_hairpin$Fg),
min(dist_multibranch$Fg), min(dist_bulge$Fg))) - 0.1
y_lim <- c(min_y, max_y)
}
plot(dist_bulge$x, Bkg$Bkg, type = "l",
col = "black", xlim = x_lim, ylim = y_lim,
main = NULL, xlab = "Nucleotide", ylab = "Density x 100",
cex.axis = 0.8, cex.lab = 1, cex.main = 1.2, lwd = 2)
arrows(dist_bulge$x, Bkg$Bkg - (error * Bkg$Bkg_se),
dist_bulge$x, Bkg$Bkg + (error * Bkg$Bkg_se),
length = 0, angle = 90, code = 3,
col = rgb(blue = 0, red = 0, green = 0, alpha = 0.5))
abline(v = 0, col = "grey", lty = 2)
lines(dist_bulge$x, dist_bulge$Fg, col = "blue", lwd = 2)
lines(dist_multibranch$x, dist_multibranch$Fg, col = "red", lwd = 2)
lines(dist_stem$x, dist_stem$Fg, col = "green", lwd = 2)
lines(dist_hairpin$x, dist_hairpin$Fg, col = "purple", lwd = 2)
lines(dist_internal$x, dist_internal$Fg, col = "orange", lwd = 2)
lines(dist_exterior$x, dist_exterior$Fg, col = "cyan", lwd = 2)
arrows(dist_bulge$x, dist_bulge$Fg - (error * dist_bulge$Fg_se),
dist_bulge$x, dist_bulge$Fg + (error * dist_bulge$Fg_se),
length = 0, angle = 90, code = 3,
col = rgb(blue = 1, red = 0, green = 0, alpha = 0.5))
arrows(dist_multibranch$x, dist_multibranch$Fg -
(error * dist_multibranch$Fg_se), dist_multibranch$x,
dist_multibranch$Fg + (error * dist_multibranch$Fg_se),
length = 0, angle = 90, code = 3,
col = rgb(blue = 0, red = 1, green = 0, alpha = 0.5))
arrows(dist_stem$x, dist_stem$Fg - (error * dist_stem$Fg_se),
dist_stem$x, dist_stem$Fg + (error * dist_stem$Fg_se),
length = 0, angle = 90, code = 3,
col = rgb(blue = 0, red = 0, green = 1, alpha = 0.5))
arrows(dist_hairpin$x, dist_hairpin$Fg - (error * dist_hairpin$Fg_se),
dist_hairpin$x, dist_hairpin$Fg + (error * dist_hairpin$Fg_se),
length = 0, angle = 90, code = 3,
col = rgb(blue = 1, red = 1, green = 0, alpha = 0.5))
arrows(dist_internal$x, dist_internal$Fg -
(error*dist_internal$Fg_se), dist_internal$x, dist_internal$Fg +
(error * dist_internal$Fg_se), length = 0, angle = 90,
code = 3, col = rgb(blue = 0, red = 1, green = 1, alpha = 0.5))
arrows(dist_exterior$x, dist_exterior$Fg -
(error * dist_exterior$Fg_se), dist_exterior$x,
dist_exterior$Fg +
(error * dist_exterior$Fg_se), length = 0, angle = 90,
code = 3, col = rgb(blue = 1, red = 0, green = 1, alpha = 0.5))
if (legend == TRUE) {
legend("bottomright",
legend = c("Background", "Bulge", "Multibranch", "Stem",
"Hairpin", "Internal", "Exterior"),
col = c("black", "blue", "red", "green", "purple",
"orange", "cyan"), lty = 1, cex = 0.8)
}
dev.off()
} else {
out_file <- paste0(out_file, ".jpeg")
coords_1 <- which(dist_stem$x == x_lim[1])
coords_2 <- which(dist_stem$x == x_lim[2])
heatmap_plot <- data.frame(stem = dist_stem$Fg[coords_1:coords_2],
hairpin = dist_hairpin$Fg[coords_1:coords_2],
bulge = dist_bulge$Fg[coords_1:coords_2],
internal = dist_internal$Fg[coords_1:coords_2],
multibranch = dist_multibranch$Fg[coords_1:coords_2],
exterior = dist_exterior$Fg[coords_1:coords_2])
heatmap_plot <- as.matrix(t(heatmap_plot))
zero <- which(dist_stem$x == 0) - coords_1
if(zero < 0){
colsep <- 0
} else{
colsep <- c(zero - 1, zero)
}
if (is.null(y_lim)) {
max_y <- max(c(max(dist_exterior$Fg), max(dist_internal$Fg),
max(dist_stem$Fg), max(dist_hairpin$Fg),
max(dist_multibranch$Fg), max(dist_bulge$Fg))) + 0.1
min_y <- min(c(min(dist_exterior$Fg), min(dist_internal$Fg),
min(dist_stem$Fg), min(dist_hairpin$Fg),
min(dist_multibranch$Fg), min(dist_bulge$Fg))) - 0.1
y_lim <- c(-max(abs(max_y), abs(min_y)),
max(abs(max_y), abs(min_y)))
}
key <- TRUE
if (legend == FALSE) {key <- FALSE}
# make plot
jpeg(out_file, height = 4, width = 8, units = "in", res = 300)
heatmap.2(heatmap_plot, Rowv = NA, Colv = NA,
col = redblue(256), dendrogram = "none",
labCol = FALSE, trace = "none", density.info = "none",
symkey = FALSE, key = key, colsep = colsep,
sepcolor = "grey", margins = c(1, 8),
breaks = seq(y_lim[1], y_lim[2], length.out = 257))
dev.off()
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.