Skip to content

Commit

Permalink
virtual4C: neater
Browse files Browse the repository at this point in the history
  • Loading branch information
Robin Van der Weide committed Oct 3, 2019
1 parent bfaff0e commit 392042b
Show file tree
Hide file tree
Showing 4 changed files with 84 additions and 47 deletions.
104 changes: 67 additions & 37 deletions R/methods_visualise.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
# Documentation -----------------------------------------------------------

#' @name visualise
#' @title Visualise discoveries
#'
Expand All @@ -24,11 +23,20 @@
#'
#' @param raw A \code{logical} of length 1: should a bare bones plot be
#' returned?
#' @param mode (PESCAn_discovery and ARA_discovery only) What result slot should
#'
#' @param bins \code{[virtual4C]} Set the number of histogram-bins
#'
#' @param bed \code{[virtual4C]} A data.frame of bed-entries to plot
#' underneath.
#'
#' @param extend_viewpoint \code{[virtual4C]} Add a set bp to both sides
#' of the viewpoint. Will make the viewpoint-box broader.
#'
#' @param mode \code{[PESCAn & ARA]} What result slot should
#' be used for visualisation? \code{"obsexp"} for the observed over expected
#' metric or \code{"signal"} for mean contacts at unshifted anchors.
#'
#' @param flipFacet (RCP_discovery only) Do you want to have RCP's of different
#' @param flipFacet \code{[RCP]} Do you want to have RCP's of different
#' regions in one plot, instead of facets? (default : \code{FALSE})
#'
#' @details The \code{"diff"} \code{metric} value creates contrast panels by
Expand Down Expand Up @@ -726,64 +734,86 @@ visualise.RCP_discovery = function(discovery, contrast = 1, metric = c("smooth",

#' @rdname visualise
#' @export
visualise.virtual4C_discovery <- function(discovery, bins = NULL, bed = NULL, bp_blackout = NULL){
visualise.virtual4C_discovery <- function(discovery, bins = NULL, bed = NULL, extend_viewpoint = NULL){
# ! someday: allow mulitple samples

data <- discovery$data
VP <- attr(discovery,"viewpoint")

if(!is.null(extend_viewpoint)){
VP[1,2] <- VP[1,2] - extend_viewpoint
VP[1,3] <- VP[1,3] + extend_viewpoint
}

if( is.null(bins) ) {
bins = nrow(data)
}

draw_blackout = F
if( is.null(bp_blackout) ) {
bp_blackout <- attr(discovery, 'resolution') * 25
draw_blackout = T
} else if (bp_blackout == F){
bp_blackout <- 0
}
VP_mid <- rowMeans(attr(discovery, 'viewpoint')[1,2:3])

blackout_up <- attr(discovery, 'viewpoint')[1,2] - (bp_blackout/2)
blackout_down <- attr(discovery, 'viewpoint')[1,3] + (bp_blackout/2)
blackout_up <- VP[1,2]
blackout_down <- VP[1,3]
data_blackout <- data[(data$mid >= blackout_up & data$mid <= blackout_down)]

if( !is.null(bed)) {
bed = bed[bed[,1] == attr(discovery, 'viewpoint')[1,1],2:3]
}

data <- data[!(data$mid > blackout_up & data$mid < blackout_down)]
data <- data[!(data$mid >= blackout_up & data$mid <= blackout_down)]

breaks <- seq(min(data$mid), max(data$mid), length.out = bins)
smooth <- data[, mean(signal),by = findInterval(data$mid, breaks)]
smooth$mid = breaks[unlist(smooth[,1])] + unique(diff(breaks)/2)
breaks <- seq(min(data$mid), max(data$mid), length.out = bins)
bin_size <- median(diff(breaks))
smooth <- data[, mean(signal),by = findInterval(data$mid, breaks)]
smooth$mid = breaks[unlist(smooth[,1])] +(bin_size/2)
smooth[,1] = NULL
colnames(smooth) = c("signal","mid")

p = ggplot2::ggplot(data, ggplot2::aes(x= mid/1e6, y = signal)) +
ggplot2::geom_col(data = smooth, fill = 'black') +
ggplot2::coord_cartesian(expand = F) +
ggplot2::geom_col(data = smooth, fill = 'black', width = bin_size/1e6) +
ggplot2::theme_classic() +
ggplot2::labs(x = attr(discovery, 'viewpoint')[1,1])

if( draw_blackout ){
p = p + ggplot2::annotate('rect',
fill = "#D8D8D8",
xmin = blackout_up/1e6,
xmax = blackout_down/1e6,
ymin = 0,
ymax = ceiling(max(smooth$signal)))
}
# draw_blackout ===+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ymax <- ceiling(max(smooth$signal))
data_blackout$signal <- ymax

p = p + ggplot2::annotate('rect',
fill = "#D8D8D8",
xmin = (min(data_blackout$mid)/1e6)-(bin_size/1e6),
xmax = (max(data_blackout$mid)/1e6)+(bin_size/1e6),
ymin = 0,
ymax = max(data_blackout$signal) )+
ggplot2::annotate(geom = 'text',
vjust = 1,
x = rowMeans(VP[,2:3])/1e6,
y = ymax*0.9,
label = '\u2693')


if( !is.null(bed)){
p = p + ggplot2::annotate('rect',
fill = "black",
xmin = bed[,1]/1e6,
xmax = bed[,2]/1e6,
ymin = -ceiling(max(smooth$signal))/100,
ymax = 0)
}
p
bed_track_size <- ceiling(max(smooth$signal))/50
bed_track_padding <- ceiling(max(smooth$signal))/200

if( !is.null(bed) ){
bed <- bed/1e6

p = p+ ggplot2::annotate(geom = 'rect',
xmin = bed[,1],
xmax = bed[,2],
ymin = -bed_track_size,
ymax = -bed_track_padding,
fill = 'black') +
ggplot2::coord_cartesian(expand = F,
ylim = c(-1*(bed_track_size+
bed_track_padding+
bed_track_padding),
ymax))
} else {
p = p + ggplot2::coord_cartesian(expand = F,
ylim = c(0, ymax))
}

p <- p + ggplot2::theme(axis.line = ggplot2::element_line(colour = 'black'),
axis.text = ggplot2::element_text(colour = 'black'))
suppressWarnings(p)

}

Expand Down
11 changes: 5 additions & 6 deletions R/virtual_4C.R
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
#' Extract matrices around a defined region a.k.a. the viewpoint
#'
#' @param exp The Hi-C experiment object of a sample: produced by construct.experiment().
#' @param viewpoint The viewpoint
#' @param viewpoint The viewpoint: will take the middle Hi-C bin if it spans mulitple bins.
#' @param xlim A vector of two with the flanking basepairs up- and downstream
#' of the viewpoint, resp.
#' @return A virtual4C_discovery object.
Expand All @@ -21,10 +21,9 @@ virtual_4C <- function(exp, viewpoint, xlim = NULL){
downstream_signal <- signal[V1 %in% vp_idx][,2:3]
} else {
# run for a region =========================================================
flank_downstream <- floor((viewpoint[,2] - xlim[1])/attr(exp, 'resolution'))
flank_upstream <- floor((xlim[2] - viewpoint[,3])/attr(exp, 'resolution'))

range_idx <- unlist(vp_idx - flank_upstream):unlist(vp_idx + flank_upstream)
flank <- floor(xlim/attr(exp, 'resolution'))

range_idx <- unlist(vp_idx - flank[1]):unlist(vp_idx + flank[2])

upstream_signal <- exp$MAT[list(range_idx,vp_idx), nomatch = 0][,c(1,3)]
downstream_signal <- exp$MAT[list(vp_idx, range_idx), nomatch = 0][,2:3]
Expand All @@ -46,7 +45,7 @@ virtual_4C <- function(exp, viewpoint, xlim = NULL){


# output ---------------------------------------------------------------------
signal <- signal[,c(3,6,2)]
signal <- unique(signal[,c(3,6,2)])
colnames(signal) <- c('chromosome','mid','signal')
signal <- list(data = signal)

Expand Down
2 changes: 1 addition & 1 deletion man/virtual_4C.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

14 changes: 11 additions & 3 deletions man/visualise.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

0 comments on commit 392042b

Please sign in to comment.