R/methods_plot.R

In robinweide/GENOVA: GENOVA: expore the Hi-C

# Documentation -----------------------------------------------------------

#' @name plot_discovery
#' @title Plot discoveries
#'
#' @description The base R plotting equivalent to
#'   \code{\link[GENOVA]{visualise}} for discovery objects.
#'
#' @param x A \code{\link[GENOVA]{discovery}} object.
#' @param contrast An \code{integer} or \code{character} matching an experiment
#'   (name) of length 1, specifying which sample should be used to make a
#'   contrast with all other samples. Alternatively, set to \code{NULL} to not
#'   plot contrast panels.
#' @param colour_fun,colour_fun_contrast A \code{function} that takes an integer
#'   as input and returns colours as output, to colour raster-like data or the
#'   constrast panels of raster-like data. Typically a function generated by
#'   \code{\link[grDevices]{colorRampPalette}}.
#' @param colour_lim,colour_lim_contrast A \code{numeric} of length 2 indicating
#'   the minimum and maximum values for which colours should be plotted.
#' @param chr A \code{character} vector of length 1 with the name of a
#'   chromosome.
#' @param start An \code{integer} of length 1 with the start basepair of the
#'   region to be plotted.
#' @param end An \code{integer} of length 1 with the ending basepair of the
#'   region to be plotted.
#' @param censor_vp \code{[virtual_4C]} A \code{logical} of length 1 deciding
#'   wether the values underneath the viewpoint should be excluded.
#' @param censor_contrast \code{[IIT]} A \code{logical} of length 1 deciding
#'   wether the contrasting experiment itself should be censored (\code{TRUE})
#'   or included (\code{FALSE}).
#' @param geom A \code{character} of length 1 indicating what geometry should be
#'   plotted. Either \code{"boxplot"} or \code{"point"}.
#' @param minimalist \code{[domainogram]} A \code{logical} of length 1. If
#'   \code{TRUE}, plot only the first experiment in the discovery and don't plot
#'   any axis, titles or legends.
#' @param metric \code{[RCP]} Currently not in use.
#' @param ... Not currently used for discovery plots.
#'
#' @return Nothing, but outputs a plot to the graphics device.
#'
#' @examples
#' \dontrun{
#' plot(discovery)
#' }
NULL

# ARMLA -------------------------------------------------------------------

#' @rdname plot_discovery
#' @export
plot.APA_discovery <- function(
  x, 
  contrast = 1,
  colour_fun = NULL,
  colour_fun_contrast = NULL,
  colour_lim = NULL,
  colour_lim_contrast = NULL,
  ...
) {
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  
  par(mar = c(1, 1, 1, 1))
  par(oma = c(4, 4, 1, 3))
  
  obj <- x$signal
  
  titles <- dimnames(obj)[[3]]
  
  xbreaks <- rev(as.numeric(dimnames(obj)[[1]])) / 1000
  ybreaks <- as.numeric(dimnames(obj)[[2]]) / 1000
  
  n_samples <- dim(obj)[3]
  n_rows <- length(contrast) + 1
  
  layout_mat <- matrix(seq_len((n_samples + 1) * n_rows), 
                       nrow = n_rows, byrow = 2)
  layout(layout_mat, widths = c(rep.int(1, n_samples), 0.2),
         respect = TRUE)
  
  if (!is.function(colour_fun)) {
    cols <- .choose_palette(colour_fun)
    colour_fun <- colorRampPalette(cols)
  }

  if (is.null(colour_lim)) {
    colour_lim <- range(obj)
  }
  replace <- which(is.na(colour_lim))
  colour_lim[replace] <- range(obj)[replace]
  
  for (i in seq_len(dim(obj)[3])) {
    m <- obj[, , i]
    m <- pmax(m, colour_lim[1])
    m <- pmin(m, colour_lim[2])
    m[cbind(rev(row(m)[T]), col(m)[T])] <- m
    image(x = xbreaks, y = ybreaks, z = m, 
          col = colour_fun(255),
          zlim = colour_lim, 
          xlab = c("Distance 3' (kb)"),
          ylab = c("Distance 5' (kb)"),
          xaxt = if(nrow(layout_mat) > 1) "n" else "s",
          yaxt = if(i == 1) "s" else "n")
    mtext(titles[[i]], side = 3, line = 0.5)
  }
  
  # A legend
  m <- t(as.matrix(seq(colour_lim[1], colour_lim[2], length.out = 255)))
  image(1, m[1, ], m, col = colour_fun(255), 
        xaxt = "n", yaxt = "n", new = FALSE)
  axis(side = 4, lwd = 0, lwd.ticks = 1, lend = 1)
  mtext("Mean Contacts", side = 4, line = 2.5)

  if (length(contrast) > 0) {
    diff <- obj
    diff[] <- diff - rep(as.vector(obj[, , contrast]), dim(diff)[3])
    
    if (is.null(colour_fun_contrast)) {
      cols <- bezier_corrected_divergent
      colour_fun_contrast <- colorRampPalette(cols)
    }
    
    if (is.null(colour_lim_contrast)) {
      colour_lim_contrast <- max(abs(range(diff))) * c(-1, 1)
    }
    replace <- which(is.na(colour_lim_contrast))
    colour_lim_contrast[replace] <- range(diff)[replace]
    
    for (i in seq_len(dim(diff)[3])) {
      m <- diff[, , i]
      m[cbind(rev(row(m)[T]), col(m)[T])] <- m
      image(x = xbreaks, y = ybreaks, z = m, 
            col = colour_fun_contrast(255),
            zlim = colour_lim_contrast, 
            xlab = c("Distance 3' (kb)"),
            ylab = c("Distance 5' (kb)"),
            yaxt = if(i == 1) "s" else "n")
    }
    
    # A legend
    m <- t(as.matrix(seq(colour_lim_contrast[1], colour_lim_contrast[2], 
                         length.out = 255)))
    image(1, m[1, ], m, col = colour_fun_contrast(255), 
          xaxt = "n", yaxt = "n", new = FALSE)
    axis(side = 4, lwd = 0, lwd.ticks = 1, lend = 1)
    
    mtext("Difference", side = 4, line = 2.5)
  }

  mtext("Distance 3' (kb)", side = 1, outer = TRUE, line = 2)
  mtext("Distance 5' (kb)", side = 2, outer = TRUE, line = 2)
}

#' @rdname plot_discovery
#' @export
plot.CSCAn_discovery <- function(
  x,
  colour_fun = NULL,
  colour_lim = NULL,
  ...
) {
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  par(mar = c(1, 1, 1, 1))
  par(oma = c(4, 4, 1, 3))
  
  obj <- x$obsexp
  
  titles_1 = dimnames(obj)[[3]]
  titles_2 = dimnames(obj)[[4]]
  
  xbreaks <- rev(as.numeric(dimnames(obj)[[1]])) / 1000
  ybreaks <- as.numeric(dimnames(obj)[[2]]) / 1000
  
  n_samples <- dim(obj)[4]
  n_rows <- dim(obj)[3] - 1
  n_cols <- n_rows * n_samples
  
  layout_mat <- matrix(seq_len((n_cols + 1) * n_rows),
                       nrow = n_rows, byrow = 2)
  layout_mat[, ncol(layout_mat)] <- layout_mat[1, ncol(layout_mat)]
  layout(layout_mat, widths = c(rep.int(1, n_cols), 0.2),
         respect = TRUE)
  
  if (is.null(colour_fun)) {
    cols <- bezier_corrected_divergent
    colour_fun <- colorRampPalette(cols)
  }
  
  if (is.null(colour_lim)) {
    colour_lim <- c(-1, 1) * max(abs(range(obj) - 1)) + 1
  }
  replace <- which(is.na(colour_lim))
  colour_lim[replace] <- range(obj)[replace]
  
  grps <- strsplit(as.character(dimnames(obj)[[3]]), "-")
  left  <- vapply(grps, `[`, character(1), 1)
  right <- vapply(grps, `[`, character(1), 2)
  com <- expand.grid(left = unique(left), right = unique(right),
                     exp = seq_len(n_samples))
  com$col <- as.numeric(interaction(com$right, com$exp))
  com$row <- as.numeric(com$left)
  com$grp <- match(with(com, paste0(left, "-", right)), dimnames(obj)[[3]])
  com <- com[order(com$row, com$col),]

  col <- 0
  for (i in seq_len(nrow(com))) {
    entry <- com[i, ]
    if (col > entry[["col"]]) {
      m <- t(as.matrix(seq(colour_lim[1], colour_lim[2], length.out = 255)))
      image(1, m[1, ], m, col = colour_fun(255), 
            xaxt = "n", yaxt = "n", new = FALSE)
      axis(side = 4, lwd = 0, lwd.ticks = 1, lend = 1)
      mtext("Observed / Expected", side = 4, line = 2.5)
      # plot.new()
      col <- entry[["col"]]
    }
    row <- entry[["row"]]
    col <- entry[["col"]]
    j <- entry[["grp"]]
    k <- entry[["exp"]]
    xaxt <- if (entry[["row"]] == max(com$row)) "s" else "n"
    yaxt <- if (entry[["col"]] == 1) "s" else "n"
    if (is.na(j)) {
      image(x = xbreaks, y = ybreaks, 
            z = matrix(0, length(xbreaks), length(ybreaks)), 
            col = "white", xaxt = xaxt, yaxt = yaxt)
      if (col == 1) {
        mtext(entry[["left"]], side = 2, line = 2)
      }
      next()
    }
    m <- obj[, , j, k]
    m <- pmax(m, colour_lim[1])
    m <- pmin(m, colour_lim[2])
    m[cbind(rev(row(m)[T]), col(m)[T])] <- m
    image(x = xbreaks, y = ybreaks, z = m,
          col = colour_fun(255),
          zlim = colour_lim,
          xlab = c("Distance 3' (kb)"),
          ylab = c("Distance 5' (kb)"),
          xaxt = xaxt, yaxt = yaxt)
    if (row == 1) {
      mtext(paste0(titles_2[entry[["exp"]]], "\n",
                   entry[["right"]]), side = 3, line = 0.5)
    }
    if (col == 1) {
      mtext(entry[["left"]], side = 2, line = 2)
    }
  }
  
  mtext("Distance (kb)", side = 1, outer = TRUE, line = 2)
  mtext("Distance (kb)", side = 2, outer = TRUE, line = 2)
}


#' @rdname plot_discovery
#' @export
plot.PESCAn_discovery <- function(
  x, 
  contrast = 1,
  colour_fun = NULL,
  colour_fun_contrast = NULL,
  colour_lim = NULL,
  colour_lim_contrast = NULL,
  ...
) {
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  par(mar = c(1, 1, 1, 1))
  par(oma = c(4, 4, 1, 3))
  
  obj <- x$obsexp
  
  titles <- dimnames(obj)[[3]]
  
  xbreaks <- rev(as.numeric(dimnames(obj)[[1]])) / 1000
  ybreaks <- as.numeric(dimnames(obj)[[2]]) / 1000
  
  n_samples <- dim(obj)[3]
  n_rows <- length(contrast) + 1
  
  layout_mat <- matrix(seq_len((n_samples + 1) * n_rows), 
                       nrow = n_rows, byrow = 2)
  layout(layout_mat, widths = c(rep.int(1, n_samples), 0.2),
         respect = TRUE)
  
  if (is.null(colour_fun)) {
    cols <- bezier_corrected_divergent
    colour_fun <- colorRampPalette(cols)
  }
  
  if (is.null(colour_lim)) {
    colour_lim <- c(-1, 1) * max(abs(range(obj) - 1)) + 1
  }
  replace <- which(is.na(colour_lim))
  colour_lim[replace] <- range(obj)[replace]
  
  for (i in seq_len(n_samples)) {
    m <- obj[, , i]
    m <- pmax(m, colour_lim[1])
    m <- pmin(m, colour_lim[2])
    m[cbind(rev(row(m)[T]), col(m)[T])] <- m
    image(x = xbreaks, y = ybreaks, z = m, 
          col = colour_fun(255),
          zlim = colour_lim, 
          xlab = c("Distance 3' (kb)"),
          ylab = c("Distance 5' (kb)"),
          xaxt = if(nrow(layout_mat) > 1) "n" else "s",
          yaxt = if(i == 1) "s" else "n")
    mtext(titles[[i]], side = 3, line = 0.5)
  }
  
  # A legend
  m <- t(as.matrix(seq(colour_lim[1], colour_lim[2], length.out = 255)))
  image(1, m[1, ], m, col = colour_fun(255), 
        xaxt = "n", yaxt = "n", new = FALSE)
  axis(side = 4, lwd = 0, lwd.ticks = 1, lend = 1)
  mtext("Observed / Expected", side = 4, line = 2.5)
  
  if (length(contrast) > 0) {
    diff <- obj
    diff[] <- diff - rep(as.vector(obj[, , contrast]), dim(diff)[3])
    
    if (is.null(colour_fun_contrast)) {
      cols <- bezier_corrected_greenPink
      colour_fun_contrast <- colorRampPalette(cols)
    }
    
    if (is.null(colour_lim_contrast)) {
      colour_lim_contrast <- max(abs(range(diff))) * c(-1, 1)
    }
    replace <- which(is.na(colour_lim_contrast))
    colour_lim_contrast[replace] <- range(diff)[replace]
    
    for (i in seq_len(dim(diff)[3])) {
      m <- diff[, , i]
      m[cbind(rev(row(m)[T]), col(m)[T])] <- m
      image(x = xbreaks, y = ybreaks, z = m, 
            col = colour_fun_contrast(255),
            zlim = colour_lim_contrast, 
            xlab = c("Distance 3' (kb)"),
            ylab = c("Distance 5' (kb)"),
            yaxt = if(i == 1) "s" else "n")
    }
    
    # A legend
    m <- t(as.matrix(seq(colour_lim_contrast[1], colour_lim_contrast[2], 
                         length.out = 255)))
    image(1, m[1, ], m, col = colour_fun_contrast(255), 
          xaxt = "n", yaxt = "n", new = FALSE)
    axis(side = 4, lwd = 0, lwd.ticks = 1, lend = 1)
    
    mtext("Difference", side = 4, line = 2.5)
  }
  
  mtext("Distance 3' (kb)", side = 1, outer = TRUE, line = 2)
  mtext("Distance 5' (kb)", side = 2, outer = TRUE, line = 2)
}

#' @rdname plot_discovery
#' @export
plot.ATA_discovery <- function(
  x, 
  contrast = 1,
  colour_fun = NULL,
  colour_fun_contrast = NULL,
  colour_lim = NULL,
  colour_lim_contrast = NULL,
  ...
) {
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  par(mar = c(1, 1, 1, 1))
  par(oma = c(4, 4, 1, 3))
  
  obj <- x$signal
  padding <- attr(x, "padding")
  
  titles <- dimnames(obj)[[3]]
  
  xbreaks <- seq(0, 1, length.out = dim(obj)[1])
  ybreaks <- seq(0, 1, length.out = dim(obj)[2])
  
  n_samples <- dim(obj)[3]
  n_rows <- length(contrast) + 1
  
  layout_mat <- matrix(seq_len((n_samples + 1) * n_rows), 
                       nrow = n_rows, byrow = 2)
  layout(layout_mat, widths = c(rep.int(1, n_samples), 0.2),
         respect = TRUE)
  
  if (!is.function(colour_fun)) {
    cols <- .choose_palette(colour_fun)
    colour_fun <- colorRampPalette(cols)
  }
  
  if (is.null(colour_lim)) {
    # colour_lim <- range(obj)
    colour_lim <- c(min(obj), quantile(obj, 0.95))
  }
  replace <- which(is.na(colour_lim))
  colour_lim[replace] <- range(obj)[replace]
  
  for (i in seq_len(dim(obj)[3])) {
    m <- obj[, , i]
    m <- pmax(m, colour_lim[1])
    m <- pmin(m, colour_lim[2])
    m[cbind(rev(row(m)[T]), col(m)[T])] <- m
    image(x = xbreaks, y = ybreaks, z = m, 
          col = colour_fun(255),
          zlim = colour_lim, 
          xlab = c("Distance 3' (kb)"),
          ylab = c("Distance 5' (kb)"),
          xaxt = "n",
          yaxt = "n")
    mtext(titles[[i]], side = 3, line = 0.5)
    if (i == 1) {
      axis(side = 2, lwd = 0, lwd.ticks = 1, lend = 1,
           at = cumsum(c(padding - 0.5, 1)) / (2 * padding),
           labels = c("3'", "5'"))
    }
    if (nrow(layout_mat) == 1) {
      axis(side = 1, lwd = 0, lwd.ticks = 1, lend = 1,
           at = cumsum(c(padding - 0.5, 1)) / (2 * padding),
           labels = c("5'", "3'"))
    }
  }
  
  # A legend
  m <- t(as.matrix(seq(colour_lim[1], colour_lim[2], length.out = 255)))
  image(1, m[1, ], m, col = colour_fun(255), 
        xaxt = "n", yaxt = "n", new = FALSE)
  axis(side = 4, lwd = 0, lwd.ticks = 1, lend = 1)
  mtext("Mean Contacts", side = 4, line = 2.5)
  
  if (length(contrast) > 0) {
    diff <- obj
    diff[] <- diff - rep(as.vector(obj[, , contrast]), dim(diff)[3])
    
    if (is.null(colour_fun_contrast)) {
      cols <- bezier_corrected_divergent
      colour_fun_contrast <- colorRampPalette(cols)
    }
    
    if (is.null(colour_lim_contrast)) {
      colour_lim_contrast <- max(abs(range(diff))) * c(-1, 1)
    }
    replace <- which(is.na(colour_lim_contrast))
    colour_lim_contrast[replace] <- range(diff)[replace]
    
    for (i in seq_len(dim(diff)[3])) {
      m <- diff[, , i]
      m[cbind(rev(row(m)[T]), col(m)[T])] <- m
      image(x = xbreaks, y = ybreaks, z = m, 
            col = colour_fun_contrast(255),
            zlim = colour_lim_contrast, 
            xlab = c("Distance 3' (kb)"),
            ylab = c("Distance 5' (kb)"),
            yaxt = "n", xaxt = "n")
      axis(side = 1, lwd = 0, lwd.ticks = 1, lend = 1,
           at = cumsum(c(padding - 0.5, 1)) / (2 * padding),
           labels = c("5'", "3'"))
      if (i == 1) {
        axis(side = 2, lwd = 0, lwd.ticks = 1, lend = 1,
             at = cumsum(c(padding - 0.5, 1)) / (2 * padding),
             labels = c("3'", "5'"))
      }
    }
    
    # A legend
    m <- t(as.matrix(seq(colour_lim_contrast[1], colour_lim_contrast[2], 
                         length.out = 255)))
    image(1, m[1, ], m, col = colour_fun_contrast(255), 
          xaxt = "n", yaxt = "n", new = FALSE)
    axis(side = 4, lwd = 0, lwd.ticks = 1, lend = 1)
    
    mtext("Difference", side = 4, line = 2.5)
  }
  
  mtext("Borders", side = 1, outer = TRUE, line = 2)
  mtext("Borders", side = 2, outer = TRUE, line = 2)
}

#' @rdname plot_discovery
#' @export
plot.ARA_discovery <- function(
  x, 
  contrast = 1,
  colour_fun = NULL,
  colour_fun_contrast = NULL,
  colour_lim = NULL,
  colour_lim_contrast = NULL,
  ...
) {
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  par(mar = c(1, 1, 1, 1))
  par(oma = c(4, 4, 1, 3))
  
  obj <- x$obsexp
  
  titles <- dimnames(obj)[[3]]
  
  xbreaks <- rev(as.numeric(dimnames(obj)[[1]])) / 1000
  ybreaks <- as.numeric(dimnames(obj)[[2]]) / 1000
  
  n_samples <- dim(obj)[3]
  n_rows <- length(contrast) + 1
  
  layout_mat <- matrix(seq_len((n_samples + 1) * n_rows), 
                       nrow = n_rows, byrow = 2)
  layout(layout_mat, widths = c(rep.int(1, n_samples), 0.2),
         respect = TRUE)
  
  if (is.null(colour_fun)) {
    cols <- bezier_corrected_divergent
    colour_fun <- colorRampPalette(cols)
  }
  
  if (is.null(colour_lim)) {
    colour_lim <- c(-1, 1) * max(abs(range(obj) - 1)) + 1
    # colour_lim <- range(obj)
  }
  replace <- which(is.na(colour_lim))
  colour_lim[replace] <- range(obj)[replace]
  
  for (i in seq_len(dim(obj)[3])) {
    m <- obj[, , i]
    m <- pmax(m, colour_lim[1])
    m <- pmin(m, colour_lim[2])
    m[cbind(row(m)[T], rev(col(m)[T]))] <- m
    image(x = xbreaks, y = ybreaks, z = m, 
          col = colour_fun(255),
          zlim = colour_lim, 
          xlab = c("Distance 3' (kb)"),
          ylab = c("Distance 5' (kb)"),
          xaxt = if(nrow(layout_mat) > 1) "n" else "s",
          yaxt = if(i == 1) "s" else "n")
    mtext(titles[[i]], side = 3, line = 0.5)
  }
  
  # A legend
  m <- t(as.matrix(seq(colour_lim[1], colour_lim[2], length.out = 255)))
  image(1, m[1, ], m, col = colour_fun(255), 
        xaxt = "n", yaxt = "n", new = FALSE)
  axis(side = 4, lwd = 0, lwd.ticks = 1, lend = 1)
  mtext("Observed / Expected", side = 4, line = 2.5)
  
  if (length(contrast) > 0) {
    diff <- obj
    diff[] <- diff - rep(as.vector(obj[, , contrast]), dim(diff)[3])
    
    if (is.null(colour_fun_contrast)) {
      cols <- bezier_corrected_greenPink
      colour_fun_contrast <- colorRampPalette(cols)
    }
    
    if (is.null(colour_lim_contrast)) {
      colour_lim_contrast <- max(abs(range(diff))) * c(-1, 1)
    }
    replace <- which(is.na(colour_lim_contrast))
    colour_lim_contrast[replace] <- range(diff)[replace]
    
    for (i in seq_len(dim(diff)[3])) {
      m <- diff[, , i]
      m[cbind(row(m)[T], rev(col(m)[T]))] <- m
      image(x = xbreaks, y = ybreaks, z = m, 
            col = colour_fun_contrast(255),
            zlim = colour_lim_contrast, 
            xlab = c("Distance 3' (kb)"),
            ylab = c("Distance 5' (kb)"),
            yaxt = if(i == 1) "s" else "n")
    }
    
    # A legend
    m <- t(as.matrix(seq(colour_lim_contrast[1], colour_lim_contrast[2], 
                         length.out = 255)))
    image(1, m[1, ], m, col = colour_fun_contrast(255), 
          xaxt = "n", yaxt = "n", new = FALSE)
    axis(side = 4, lwd = 0, lwd.ticks = 1, lend = 1)
    
    mtext("Difference", side = 4, line = 2.5)
  }
  
  mtext("Distance (kb)", side = 1, outer = TRUE, line = 2)
  mtext("Distance (kb)", side = 2, outer = TRUE, line = 2)
}


# Genome-wide scores ------------------------------------------------------

#' @rdname plot_discovery
#' @export
plot.CS_discovery <- function(x, chr = "chr1", start = NULL, end = NULL,
                              contrast = NULL, ...) {
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  start <- if (is.null(start)) -Inf else start
  end <- if (is.null(end)) Inf else end
  loc <- standardise_location(chr, start, end, singular = TRUE)
  df <- as.data.table(x$compart_scores)
  ii <- which(df[["chrom"]] == loc$chrom & df[["start"]] >= loc$start & 
                df[["end"]] <= loc$end)
  df <- df[ii,]
  expnames <- colnames(df)[5:ncol(df)]
  dat <- as.matrix(df[, 5:ncol(df)])
  df[, "mid" := (start + end) / 2]
  
  cols <- attr(x, "colours")
  
 
  ylim <- max(abs(quantile(dat, 0.999, na.rm = TRUE))) * c(-1, 1)
  # ylim <- max(abs(range(dat, na.rm = TRUE))) * c(-1, 1)
  plot.new()
  leg <- legend(par("usr")[2] + 5, y = 0.5, expnames, col = cols, lwd = 1, 
                title = "Sample", xpd = TRUE, bty = "n", plot = FALSE)
  leg <- leg$rect$w
  
  par(mar = c(5, 4, 4, 6) + 0.1)

  plot(df$mid / 1e6, df[[expnames[[1]]]], type = 'l', col = cols[1],
       ylim = ylim, ylab = "Compartment Score", 
       xlab = paste0("Location ", chr, " (Mb)"), new = FALSE)
  for(i in tail(seq_len(ncol(dat)), -1)) {
    lines(df$mid / 1e6, dat[, i], col = cols[i])
  }
  
  legend(par("usr")[2] + 5, y = 0.5, expnames, col = cols, lwd = 1, 
         title = "Sample", xpd = TRUE, bty = "n")
}

#' @rdname plot_discovery
#' @export
plot.IS_discovery <- function(x, chr = "chr1", start = NULL, end = NULL,
                              contrast = NULL, ...) {
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  start <- if (is.null(start)) -Inf else start
  end <- if (is.null(end)) Inf else end
  loc <- standardise_location(chr, start, end, singular = TRUE)
  df <- as.data.table(x$insula_score)
  ii <- which(df[["chrom"]] == loc$chrom & df[["start"]] >= loc$start & 
                df[["end"]] <= loc$end)
  df <- df[ii,]
  expnames <- colnames(df)[5:ncol(df)]
  dat <- as.matrix(df[, 5:ncol(df)])
  df[, "mid" := (start + end) / 2]
  
  cols <- attr(x, "colours")
  
  
  ylim <- max(abs(quantile(dat, 0.999, na.rm = TRUE))) * c(-1, 1)
  # ylim <- max(abs(range(dat, na.rm = TRUE))) * c(-1, 1)
  par(mar = c(5, 4, 4, 6) + 0.1)
  plot.new()
  leg <- legend(par("usr")[2] + 5, y = 0.5, expnames, col = cols, lwd = 1, 
                title = "Sample", xpd = TRUE, bty = "n", plot = FALSE)
  leg <- leg$rect$w
  
  
  
  plot(df$mid / 1e6, df[[expnames[[1]]]], type = 'l', col = cols[1],
       ylim = ylim, ylab = "Insulation Score", 
       xlab = paste0("Location ", chr, " (Mb)"), new = FALSE)
  for(i in tail(seq_len(ncol(dat)), -1)) {
    lines(df$mid / 1e6, dat[, i], col = cols[i])
  }
  
  legend(par("usr")[2] + 5, y = 0.5, expnames, col = cols, lwd = 1, 
         title = "Sample", xpd = TRUE, bty = "n")
}

#' @rdname plot_discovery
#' @export
plot.DI_discovery <- function(x, chr = "chr1", start = NULL, end = NULL,
                              contrast = NULL, ...) {
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  start <- if (is.null(start)) -Inf else start
  end <- if (is.null(end)) Inf else end
  loc <- standardise_location(chr, start, end, singular = TRUE)
  df <- x$DI
  ii <- which(df[["chrom"]] == loc$chrom & df[["start"]] >= loc$start & 
                df[["end"]] <= loc$end)
  df <- df[ii,]
  expnames <- tail(colnames(df), -4)
  # df <- dcast(df, chrom + start + end + bin ~ experiment, value.var = "DI")
  dat <- as.matrix(df[, 5:ncol(df)])
  setDT(df)
  df[, "mid" := (start + end) / 2]
  
  cols <- attr(x, "colours")
  
  
  ylim <- max(abs(quantile(dat, 0.999, na.rm = TRUE))) * c(-1, 1)
  # ylim <- max(abs(range(dat, na.rm = TRUE))) * c(-1, 1)
  par(mar = c(5, 4, 4, 6) + 0.1)
  plot.new()
  leg <- legend(par("usr")[2] + 5, y = 0.5, expnames, col = cols, lwd = 1, 
                title = "Sample", xpd = TRUE, bty = "n", plot = FALSE)
  leg <- leg$rect$w
  
  
  
  plot(df$mid / 1e6, df$test1, type = 'l', col = cols[1],
       ylim = ylim, ylab = "Directionality Index", 
       xlab = paste0("Location ", chr, " (Mb)"), new = FALSE)
  for(i in tail(seq_len(ncol(dat)), -1)) {
    lines(df$mid / 1e6, dat[, i], col = cols[i])
  }
  
  legend(par("usr")[2] + 5, y = 0.5, expnames, col = cols, lwd = 1, 
         title = "Sample", xpd = TRUE, bty = "n")
}


# Miscellaneous -----------------------------------------------------------

#' @rdname plot_discovery
#' @export
plot.saddle_discovery <- function(
  x, 
  contrast = 1,
  colour_fun = NULL,
  colour_fun_contrast = NULL,
  colour_lim = NULL,
  colour_lim_contrast = NULL,
  ...
) {
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  par(mar = c(1, 1, 1, 1))
  par(oma = c(4, 4, 1, 3))
  
  df <- x$saddle
  df <- df[!is.na(mean) & !is.na(q1),]
  df <- df[, mean(mean), by = list(exp, q1, q2)]
  setnames(df, 4, "obsexp")
  
  obj <- array(0, dim = c(max(df$q1), 
                          max(df$q2), 
                          length(unique(df$exp))))
  dimnames(obj)[[3]] <- unique(df$exp)
  dimnames(obj)[[1]] <- dimnames(obj)[[2]] <- seq(1, max(df$q1))
  obj[cbind(df$q1, df$q2, df$exp)] <- df$obsexp
  obj[cbind(df$q2, df$q1, df$exp)] <- df$obsexp
  
  titles <- df[, unique(exp)]
  
  xbreaks <- seq(0, 1, length.out = dim(obj)[1])
  ybreaks <- seq(0, 1, length.out = dim(obj)[2])
  
  n_samples <- dim(obj)[3]
  n_rows <- length(contrast) + 1
  
  layout_mat <- matrix(seq_len((n_samples + 1) * n_rows), 
                       nrow = n_rows, byrow = 2)
  layout(layout_mat, widths = c(rep.int(1, n_samples), 0.2),
         respect = TRUE)
  
  if (is.null(colour_fun)) {
    cols <- bezier_corrected_divergent
    colour_fun <- colorRampPalette(cols)
  }
  
  if (is.null(colour_lim)) {
    colour_lim <- c(-1, 1) * max(abs(range(obj) - 1)) + 1
    # colour_lim <- range(obj)
  }
  replace <- which(is.na(colour_lim))
  colour_lim[replace] <- range(obj)[replace]
  
  for (i in seq_len(dim(obj)[3])) {
    m <- obj[, , i]
    m <- pmax(m, colour_lim[1])
    m <- pmin(m, colour_lim[2])
    # m[cbind(rev(row(m)[T]), col(m)[T])] <- m
    image(x = xbreaks, y = ybreaks, z = m, 
          col = colour_fun(255),
          zlim = colour_lim, 
          xlab = c("Distance 3' (kb)"),
          ylab = c("Distance 5' (kb)"),
          xaxt = if(nrow(layout_mat) > 1) "n" else "s",
          yaxt = if(i == 1) "s" else "n")
    mtext(titles[[i]], side = 3, line = 0.5)
  }
  
  # A legend
  m <- t(as.matrix(seq(colour_lim[1], colour_lim[2], length.out = 255)))
  image(1, m[1, ], m, col = colour_fun(255), 
        xaxt = "n", yaxt = "n", new = FALSE)
  axis(side = 4, lwd = 0, lwd.ticks = 1, lend = 1)
  mtext("Observed / Expected", side = 4, line = 2.5)
  
  if (length(contrast) > 0) {
    diff <- obj
    diff[] <- diff - rep(as.vector(obj[, , contrast]), dim(diff)[3])
    
    if (is.null(colour_fun_contrast)) {
      cols <- bezier_corrected_greenPink
      colour_fun_contrast <- colorRampPalette(cols)
    }
    
    if (is.null(colour_lim_contrast)) {
      colour_lim_contrast <- max(abs(range(diff))) * c(-1, 1)
    }
    replace <- which(is.na(colour_lim_contrast))
    colour_lim_contrast[replace] <- range(diff)[replace]
    
    for (i in seq_len(dim(diff)[3])) {
      m <- diff[, , i]
      # m[cbind(rev(row(m)[T]), col(m)[T])] <- m
      image(x = xbreaks, y = ybreaks, z = m, 
            col = colour_fun_contrast(255),
            zlim = colour_lim_contrast, 
            xlab = c("Distance 3' (kb)"),
            ylab = c("Distance 5' (kb)"),
            yaxt = if(i == 1) "s" else "n")
    }
    
    # A legend
    m <- t(as.matrix(seq(colour_lim_contrast[1], colour_lim_contrast[2], 
                         length.out = 255)))
    image(1, m[1, ], m, col = colour_fun_contrast(255), 
          xaxt = "n", yaxt = "n", new = FALSE)
    axis(side = 4, lwd = 0, lwd.ticks = 1, lend = 1)
    
    mtext("Difference", side = 4, line = 2.5)
  }
  
  mtext("Quantile", side = 1, outer = TRUE, line = 2)
  mtext("Quantile", side = 2, outer = TRUE, line = 2)
}

#' @rdname plot_discovery
#' @export
plot.domainogram_discovery <- function(
  x, colour_fun = NULL, colour_lim = c(-1, 1), 
  minimalist = FALSE,
  ...
) {
  minimalist <- literalTRUE(minimalist)
  if (minimalist) {
    opar <- par(no.readonly = TRUE)
    on.exit(par(opar))
  }
  x <- x$scores
  if (minimalist) {
    x <- x[, 1:3]
  } else {
    par(mar = c(1, 1, 1, 1))
    par(oma = c(4, 4, 1, 3))
  }
  
  x <- as.data.table(x)
  x <- melt.data.table(x, id.vars = c("window", "position"), 
                       value.name = "insulation")
  setnames(x, 3, "experiment")
  
  mats <- x
  mats <- split(mats, mats$experiment)
  mats <- lapply(mats, function(mat) {
    casted <- dcast(mat, position ~ window, value.var = "insulation")
    rows <- casted$position
    cols <- tail(colnames(casted), -1)
    mat <- as.matrix(casted[, -1])
    dimnames(mat) <- list(rows, cols)
    mat
  })
  
  titles <- names(mats)
  
  xbreaks <- rev(as.numeric(dimnames(mats[[1]])[[1]])) / 1e6
  ybreaks <- as.numeric(dimnames(mats[[1]])[[2]])
  
  n_samples <- length(mats)
  n_rows <- 1
  
  layout_mat <- matrix(c(seq_len(n_samples), rep(n_samples + 1, n_samples)), 
                       nrow = n_samples, byrow = FALSE)
  if (!minimalist) {
    layout(layout_mat, widths = c(1, 0.1),
           respect = FALSE)
  }
  
  if (is.null(colour_fun)) {
    cols <- rev(bezier_corrected_divergent)
    colour_fun <- colorRampPalette(cols)
  }
  
  if (is.null(colour_lim)) {
    colour_lim <- range(vapply(mats, range, numeric(2)))
  }
  replace <- which(is.na(colour_lim))
  colour_lim[replace] <- range(unlist(mats))[replace]
  
  for (i in seq_len(n_samples)) {
    m <- mats[[i]]
    m <- pmax(m, colour_lim[1])
    m <- pmin(m, colour_lim[2])
    # m[cbind(rev(row(m)[T]), col(m)[T])] <- m
    if (!minimalist) {
      image(x = as.numeric(dimnames(m)[[1]]) / 1e6, 
            y = as.numeric(dimnames(m)[[2]]), 
            z = m, 
            col = colour_fun(255),
            zlim = colour_lim, 
            xlab = paste0("Location ", attr(x, "chrom"), " (Mb)"),
            ylab = "Window Size",
            xaxt = if(i == nrow(layout_mat)) "s" else "n",
            yaxt = "s")
      mtext(titles[[i]], side = 3, line = 0.5)
    } else {
      image(x = as.numeric(dimnames(m)[[1]]) / 1e6,
            y = as.numeric(dimnames(m)[[2]]),
            z = m,
            col = colour_fun(255),
            zlim = colour_lim,
            xlab = "", ylab = "",
            xaxt = "n", yaxt = "n")
    }
  }
  
  # A legend
  if (!minimalist) {
    m <- t(as.matrix(seq(colour_lim[1], colour_lim[2], length.out = 255)))
    image(1, m[1, ], m, col = colour_fun(255), 
          xaxt = "n", yaxt = "n", new = FALSE)
    axis(side = 4, lwd = 0, lwd.ticks = 1, lend = 1)
    mtext("Insulation Score", side = 4, line = 2.5)
    
    mtext(paste0("Location ", attr(x, "chrom"), " (Mb)"), 
          side = 1, outer = TRUE, line = 2)
    mtext("Window Size", side = 2, outer = TRUE, line = 2)
  }
}

#' @rdname plot_discovery
#' @export
plot.RCP_discovery <- function(x, contrast = 1, 
                               metric = c("smooth"), ...) {
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  metric <- match.arg(metric, c("smooth", "both", "lfc"))
  
  nregion <- length(unique(x$smooth$region))
  
  dats <- split(x$smooth, x$smooth$samplename)
  
  xlim = range(x$smooth$distance[x$smooth$distance > 0] / 1e6)
  ylim = range(x$smooth$P)
  
  plot(dats[[1]]$distance[-1] / 1e6, 
       dats[[1]]$P[-1], 
       log = "xy", type = 'l',
       xlab = "Distance (Mb)", ylab = "RCP", col = dats[[1]]$colour[1],
       xlim = xlim, ylim = ylim)
  
  for(i in tail(seq_along(dats), -1)) {
    lines(dats[[i]]$distance[-1] / 1e6, dats[[i]]$P[-1], col = dats[[i]]$colour[1])
  }
  
}

#' @rdname plot_discovery
#' @export
plot.virtual4C_discovery <- function(x, censor_vp = TRUE, ...) {
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  
  par(mar = c(1, 1, 1, 1))
  par(oma = c(4, 4, 1, 1))
  
  res <- attr(x, "resolution")
  vp  <- attr(x, "viewpoint")
  
  data <- as.data.table(x$data)
  data <- melt.data.table(data, id.vars = c("chromosome", "mid"))
  colnames(data) <- c("chromosome", "mid", "experiment", "signal")
  data <- data[chromosome == vp[1,1, drop = TRUE]]
  
  # fill in empty signal
  sdata <- split(data, data$experiment)
  sdata <- lapply(sdata, function(dat) {
    missing <- seq(min(dat$mid), max(dat$mid), by = res)
    missing <- missing[!(missing %in% dat$mid)]
    append <- data.table(dat[1, 1], missing, dat[1, 3], 0)
    dat <- rbind(dat, append, use.names = FALSE)
    dat[order(mid)]
  })
  data <- rbindlist(sdata)
  data <- data[is.finite(signal)]
  
  data <- data[, list(chromosome, pos = mid + c(-1, 1) * res, 
              signal, experiment, p = c("s", "e")), by = seq_len(nrow(data))]
  data <- dcast(data, chromosome + pos + p ~ experiment, value.var = "signal")
  
  expnames <- tail(colnames(data), - 3)
  
  n_samples <- length(expnames)
  n_rows <- 1
  
  layout_mat <- matrix(c(seq_len(n_samples), rep(n_samples + 1, n_samples)), 
                       nrow = n_samples, byrow = FALSE)
  layout(layout_mat, widths = c(1, 0.1),
         respect = FALSE)
  
  for (i in expnames) {
    i <- as.symbol(i)
    data[, as.character(i) := pmax(eval(i), 0, na.rm = TRUE)]
  }
  
  start <- head(data, 1)
  end <- tail(data, 1)
  start[, eval(expnames) := 0]
  end[, eval(expnames) := 0]
  
  for (i in expnames) {
    i <- as.symbol(i)
    start[, as.character(i) := 0]
    end[, as.character(i) := 0]
  }
  
  if (censor_vp) {
    if (nrow(vp) == 1) {
      data[pos >= vp[1,2] & pos <= vp[1,3], eval(expnames) := 0]
    } else {
      for (i in seq_along(expnames)) {
        this <- vp[vp$exp == expnames[i], ]
        data[pos >= this[1, 2] & pos <= this[1, 3], eval(expnames[[i]]) := 0]
      }
    }
  }
  
  data <- rbind(start, data, end)
  data$pos <- data$pos / 1e6
  
  ylim <- c(0, max(data[, eval(expnames), with = FALSE]))
  xlim <- range(data$pos)
  
  plot(data$pos, data[[expnames[1]]], type = "n", ylim = ylim, xlim = xlim,
       xaxt = if(nrow(layout_mat) > 1) "n" else "s")
  polygon(data$pos, data[[expnames[[1]]]], col = "black")
  
  mtext(expnames[[1]], side = 3, line = 0.1)
  
  if (censor_vp) {
    this <- vp[vp$exp == expnames[1], ]
    rect(this[1, 2] / 1e6, ylim[1] - ylim[2], 
         this[1, 3] / 1e6, ylim[2] * 2, 
         col = "#D8D8D8", border = NA)
    text(x = mean(c(this[1, 2], this[1, 3])) / 1e6,
         y = ylim[1] + diff(ylim) * 0.80, labels = '\u2693')
  }
  
  for (i in tail(seq_along(expnames), -1)) {
    plot(data$pos, data[[expnames[i]]], type = "n", ylim = ylim, xlim = xlim,
         xaxt = if (i == nrow(layout_mat)) "s" else "n")
    polygon(data$pos, data[[expnames[i]]], col = "black")
    
    if (censor_vp) {
      if (nrow(vp) == 1) {
        rect(vp[1,2] / 1e6, ylim[1] - ylim[2], 
             vp[1,3] / 1e6, ylim[2] * 2, 
             col = "#D8D8D8", border = NA)
        text(x = mean(c(vp[1, 2], vp[1, 3])) / 1e6,
             y = ylim[1] + diff(ylim) * 0.80, labels = '\u2693')
      } else {
        this <- which(vp$exp == expnames[i])
        rect(vp[this, 2] / 1e6, ylim[1] - ylim[2],
             vp[this, 3] / 1e6, ylim[2] * 2,
             col = "#D8D8D8", border = NA)
        text(x = mean(c(vp[this, 2], vp[this, 3])) / 1e6,
             y = ylim[1] + diff(ylim) * 0.80, labels = '\u2693')
      }

    }
    mtext(expnames[[i]], side = 3, line = 0.1)
  }
  
  mtext(paste0("Location ", vp[1, 1], " (Mb)"), 
        side = 1, line = 2, outer = TRUE)
  mtext("Signal", side = 2, line = 2, outer = TRUE)
}

#' @rdname plot_discovery
#' @export
plot.IIT_discovery <- function(x, contrast = 1, censor_contrast = TRUE, 
                               geom = c("boxplot", "point"),
                               ...) {
  opar <- par(no.readonly = TRUE)
  on.exit(par(opar))
  geom <- match.arg(geom)
  dat <- as.data.table(x$results)
  cols <- attr(x, "colours")
  
  expnames <- tail(colnames(dat), -2)
  
  if (!is.null(contrast) & length(expnames) < 2) {
    message("Cannot compute a contrast for one sample. Reverting to ",
            "plotting plain values.")
    contrast <- NULL
  } else if (!is.null(contrast)) {
    contrast <- expnames[contrast]
    
    trans <- lapply(setNames(expnames, expnames), function(i) {
      log2(dat[[i]] / dat[[contrast]])
    })
    
    for (i in expnames) {
      dat[, as.character(i) := trans[[i]]]
    }
    
    if (censor_contrast & !is.null(contrast)) {
      dat <- dat[, -..contrast]
      cols <- cols[which(expnames != contrast)]
    }
  }
  
  df <- melt.data.table(
    dat, id.vars = c("x", "y"),
    measure.vars = intersect(colnames(dat), expnames)
  )
  
  df$diff <- as.factor(df$y - df$x)
  
  udiffs <- unique(df$diff)
  ndiffs <- length(udiffs)
  nvars <- length(unique(df$variable))
  
  groups <- rep(seq_len(ndiffs), 
                each = nvars)
  xpos <- groups + seq_along(groups) - 1
  xticks <- unlist(by(xpos, groups, mean, simplify = FALSE))
  
  if (is.null(contrast)) {
    ylab <- "Contacts"
    df$value <- log10(df$value)
  } else {
    ylab <- bquote("Log"[2]~" (Experiment contacts /" ~
                     paste(.(contrast)) ~ "contacts)")
  }
  
  if (geom == "boxplot") {
    boxplot(value ~ variable + diff, df,
            at = xpos, xaxt = "n", yaxt = "n",
            col = cols, pch = 19, cex = 0.5,
            xlab = "TAD Distance", ylab = ylab)
  } else {
    xjitter <- jitter(xpos[as.numeric(interaction(df$variable, df$diff))])
    plot(xjitter, 
         df$value,
         col = scales::alpha(cols[as.numeric(df$variable)], 0.3),
         pch = 19, cex = 0.5, xaxt = "n", yaxt = "n",
         xlab = "TAD Distance", ylab = ylab)
  }
  
  axis(1, at = xticks, labels = paste0("n + ", udiffs),
       lwd = 0, lwd.ticks = 1)
  
  ybreaks <- scales::extended_breaks()(range(df$value))
  if (is.null(contrast)) {
    ylabels <- scales::math_format()(ybreaks)
    axis(2, at = ybreaks, labels = as.expression(ylabels), las = 1,
         lwd = 0, lwd.ticks = 1)
  } else {
    ylabels <- scales::format_format()(ybreaks)
    axis(2, at = ybreaks, labels = ylabels, las = 1,
         lwd = 0, lwd.ticks = 1)
  }
  
  pos <- if (mean(df$value[df$diff == 1]) > mean(range(df$value))) {
    "bottomleft"
  } else {
    "topleft"
  }
  
  legend(pos, legend = intersect(colnames(dat), expnames),
         fill = cols, bty = "n")
  
}