RCPvis: ylab for lfc

loader: colour-bug and read.hicpro.matrix transfer

NAMESPACE

@@ -15,6 +15,7 @@ S3method(visualise,ATA_discovery)
 S3method(visualise,PESCAn_discovery)
 S3method(visualise,RCP_discovery)
 S3method(visualise,default)
+S3method(visualise,virtual4C_discovery)
 export(APA)
 export(APA_old)
 export(ARA)
@@ -68,6 +69,7 @@ export(select_subset)
 export(shuffleHiC)
 export(trans.compartment.plot)
 export(unbundle)
+export(virtual_4C)
 export(visualise)
 export(visualise_APA_base)
 export(visualise_APA_ggplot)

R/load_contacts.R

@@ -199,7 +199,9 @@ load_contacts = function(signal_path,
   ##############################################################################
   ###############################################################  Contruct list
   ##############################################################################
-  colour = if(is.null(colour)){colour = 'black'}
+  if( is.null(colour) ){
+    colour = 'black'
+  }
   out = structure(list(
     # Iced HiC-matrix in three-column format (i.e. from HiC-pro)
     MAT = signal,
@@ -356,3 +358,22 @@ legacy_out = function(out){
   return(newOut)
 
 }
+
+#' Read hicpro three column matrix format.
+#'
+#' This function loads a HiC-pro file as a matrix. It assumes a three-column
+#' layout: bin1, bin2 and score. All scors are normalised to contacts per *scale_bp*
+#' total contacts.
+#'
+#' @param file Full path to file.
+#' @param scale_bp Normalising factor. Set to NULL to skip scale_bp.
+#' @return A data.table with normalised counts.
+read.hicpro.matrix <- function(file, scale_bp = 1e9) {
+  data <- data.table::fread(file)
+  data.table::setkey(data, "V1", "V2")
+  if (!is.null(scale_bp)) {
+    data$V3 <- scale_bp * data$V3 / sum(data$V3)
+  }
+  return(data)
+}
+

R/methods_visualise.R

@@ -612,20 +612,23 @@ visualise.RCP_discovery = function(discovery, contrast = 1, metric = c("smooth",
     })
     lfcDT = rbindlist(lfcDT)
 
+
+
+
     GG = NULL
     if(nregion == 1){
       GG = ggplot2::ggplot(lfcDT, ggplot2::aes(x = log10(distance), y = P ,col = samplename)) +
-        ggplot2::labs(x = 'distance (Mb)', col = 'sample') +
+        ggplot2::labs(x = 'distance (Mb)', col = 'sample', y = expression("log2(P"[sample]*"/P"[contrast]*")")) +
         ggplot2::scale_color_manual(values = smplCols)
     } else if(flipFacet){
       GG =   ggplot2::ggplot(lfcDT, ggplot2::aes(x = log10(distance), y = P ,col = region)) +
         ggplot2::facet_grid(. ~ samplename)+
-        ggplot2::labs(x = 'distance (Mb)', col = 'region') +
+        ggplot2::labs(x = 'distance (Mb)', col = 'region', y = expression("log2(P"[sample]*"/P"[contrast]*")")) +
         ggplot2::scale_color_manual(values = regCols)
     } else {
       GG =ggplot2::ggplot(lfcDT, ggplot2::aes(x = log10(distance), y = P ,col = samplename)) +
         ggplot2::facet_grid(. ~ region)+
-        ggplot2::labs(x = 'distance (Mb)', col = 'sample')+
+        ggplot2::labs(x = 'distance (Mb)', col = 'sample', y = expression("log2(P"[sample]*"/P"[contrast]*")"))+
         ggplot2::scale_color_manual(values = smplCols)
     }
     RAUW = GG +   
@@ -720,6 +723,70 @@ 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){
+  # ! someday: allow mulitple samples
+
+  data <- discovery$data
+  VP   <- attr(discovery,"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
+  } 
+
+  blackout_up   <- attr(discovery, 'viewpoint')[1,2] - (bp_blackout/2)
+  blackout_down <- attr(discovery, 'viewpoint')[1,3] + (bp_blackout/2)
+
+  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)]
+
+  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)
+  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::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))) 
+  }
+
+  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
+
+
+}
+
 # Utilities ---------------------------------------------------------------
 
 # Makes sure no errors are returned when visualise(..., raw = TRUE)

R/virtual_4C.R

@@ -1,13 +1,15 @@
-VP = HICCUPS[2498,1:3]
-xlim      <- VP[,2:3]
-xlim[1]   <- xlim[1]-5e6
-xlim[2]   <- xlim[2]+5e6
-
-virtual_4C(exp = WT_new, viewpoint = VP, xlim = xlim)
-
-
+#' Get an average virtual 4C of 2D-regions (e.g. loops).
+#'
+#' 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 xlim A vector of two with the flanking basepairs up- and downstream 
+#' of the viewpoint, resp. 
+#' @return A virtual4C_discovery object.
+#' @export
 virtual_4C <- function(exp, viewpoint, xlim = NULL){
-
+  # ! someday: allow mulitple samples
   vp_idx <- median(bed2idx(exp$IDX, viewpoint))
 
   signal <- NULL
@@ -46,77 +48,17 @@ virtual_4C <- function(exp, viewpoint, xlim = NULL){
   # output ---------------------------------------------------------------------
   signal <- signal[,c(3,6,2)]
   colnames(signal) <- c('chromosome','mid','signal')
+  signal <- list(data = signal)
 
   signal <- structure(signal, 
+                      class = "virtual4C_discovery",
                       'viewpoint' = viewpoint, 
                       'xlim' = xlim,
                       'sample' = attr(exp, 'sample'),
-                      class = "virtual4C_discovery",
+                      'resolution' = attr(exp, 'resolution'),
                       package = "GENOVA")
 
   signal
 }
 
 
-
-visualise.virtual4C_discovery(){
-
-
-}
-
-
-library(ggplot2)
-bed = HICCUPS[,1:3]
-bins = 200
-blackout_region <- 2.5e5
-data <- signal
-draw_blackout <- T
-
-if( is.null(bins) ) {
-  bins = nrow(data)
-}
-
-if( !is.null(bed)) {
-  bed = bed[bed[,1] == attr(data, 'viewpoint_chromosome'),2:3]
-}
-
-blackout_up   <- attr(data, 'viewpoint_start') - (blackout_region/2)
-blackout_down <- attr(data, 'viewpoint_end') + (blackout_region/2)
-
-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)
-smooth[,1] = NULL
-colnames(smooth) = c("signal","mid")
-
-p = ggplot(data, aes(x= mid/1e6, y = signal)) +
-  annotate('rect', 
-           fill = "black",
-           xmin = bed[,1]/1e6, 
-           xmax = bed[,2]/1e6,
-           ymin = -ceiling(max(smooth$signal))/100,
-           ymax = 0)  +
-  geom_col(data = smooth, fill = 'black') +
-  coord_cartesian(xlim = unlist(attr(signal, 'xlim'))/1e6 ,expand = F) +
-  theme_classic() 
-
-if( draw_blackout ){
-  p = p + annotate('rect',
-                   fill =  "#D8D8D8",
-                   xmin = blackout_up/1e6,
-                   xmax = blackout_down/1e6,
-                   ymin = 0,
-                   ymax = ceiling(max(smooth$signal))) 
-}
-
-p
-
-
-
-
-
-
-
-