Added multiple experiments as option for virtual 4C, see robinweide#123…

…. Also added print method.

R/methods_print.R

@@ -280,6 +280,35 @@ print.domainogram_discovery <- function(x, ...) {
   print(as.data.table(x))
 }
 
+#' @export
+#' @keywords internal
+print.virtual4C_discovery <- function(x, ...) {
+  myclass <- class(x)[[1]]
+  expnames <- unique(x$data$experiment)
+  res <- attr(x, "resolution")
+  res <- if (res %% 1e6 == 0) {
+    paste0(res / 1e6, " Mb")
+  } else if (res %% 1e3 == 0) {
+    paste0(res / 1e3, " kb")
+  } else {
+    paste0(res, " bp")
+  }
+
+  vp <- attr(x, "viewpoint")
+
+  string <- paste0("A ", attr(x, "package"), " '", myclass,
+                   "' object involving the following ",
+                   length(expnames), " experiments:\n'",
+                   paste0(expnames, collapse = "', '"), "' at a ",
+                   "resolution of ", res, ".\n")
+  string1 <- paste0("The viewpoint of this virtual 4C is located at ",
+                    vp[1, 1], ":", format(vp[1, 2], scientific = FALSE), "-",
+                    format(vp[1, 3], scientific = FALSE), ".")
+
+  cat(string)
+  cat(string1)
+}
+
 # Other classes -----------------------------------------------------------
 
 #' @export

R/methods_visualise.R

@@ -49,8 +49,6 @@
 #'   \code{-Inf} and \code{Inf} respectively.
 #'   
 #' @param title add a title
-#' @param signal_size The width/height of the signal (e.g. a value of 3 will 
-#' take the middle 3x3 matrix of the APA).
 #' @param ... further arguments passed to or from other methods.
 #'
 #' @details The \code{"diff"} \code{metric} value creates contrast panels by
@@ -1013,7 +1011,7 @@ visualise.virtual4C_discovery <- function(discovery, bins = NULL, bed = NULL, ex
 
   blackout_up   <- VP[1,2]
   blackout_down <- VP[1,3]
-  data_blackout <- data[(data$mid >= blackout_up & data$mid <= blackout_down)]
+  data_blackout <- data[mid >= blackout_up & mid <= blackout_down]
 
   if( !is.null(bed)) {
     bed = bed[bed[,1] == attr(discovery, 'viewpoint')[1,1],2:3]
@@ -1023,31 +1021,36 @@ visualise.virtual4C_discovery <- function(discovery, bins = NULL, bed = NULL, ex
 
   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   <- 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', width = bin_size/1e6) +
     ggplot2::theme_classic() +
-    ggplot2::labs(x = attr(discovery, 'viewpoint')[1,1])
+    ggplot2::labs(x = attr(discovery, 'viewpoint')[1,1]) +
+    ggplot2::facet_grid(experiment ~ .)
 
   # 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 (nrow(data_blackout) > 0) {
+    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', 
@@ -1081,7 +1084,8 @@ visualise.virtual4C_discovery <- function(discovery, bins = NULL, bed = NULL, ex
   }
 
   p <- p + ggplot2::theme(axis.line = ggplot2::element_line(colour = 'black'),
-                          axis.text = ggplot2::element_text(colour = 'black'))
+                          axis.text = ggplot2::element_text(colour = 'black'),
+                          strip.background = element_blank())
   suppressWarnings(p)
 }
 

R/virtual_4C.R

@@ -8,9 +8,17 @@
 #' of the viewpoint, resp. 
 #' @return A virtual4C_discovery object.
 #' @export
-virtual_4C <- function(exp, viewpoint, xlim = NULL){
+virtual_4C <- function(explist, viewpoint, xlim = NULL){
   # ! someday: allow mulitple samples
-  vp_idx <- median(bed2idx(exp$IDX, viewpoint))
+  explist  <- GENOVA:::check_compat_exp(explist)
+  expnames <- if (is.null(names(explist))) {
+    vapply(explist, attr, character(1L), "samplename")
+  } else {
+    names(explist)
+  }
+  IDX <- explist[[1]]$IDX
+
+  vp_idx <- median(bed2idx(IDX, viewpoint))
 
   # Restrict data.table core usage
   dt.cores <- data.table::getDTthreads()
@@ -20,38 +28,47 @@ virtual_4C <- function(exp, viewpoint, xlim = NULL){
   signal <- NULL
   if( is.null(xlim) ){
     # run genome-wide ==========================================================
-    signal <- exp$MAT[V1 == vp_idx | V2 == vp_idx]
-
-    upstream_signal   <- signal[V2 %in% vp_idx][,c(1,3)]
-    downstream_signal <- signal[V1 %in% vp_idx][,2:3]
+    signal <- lapply(seq_along(explist), function(i) {
+      explist[[i]]$MAT[V1 == vp_idx | V2 == vp_idx, 
+                       list(V1, V2, V3, exp = i)]
+    })
+    signal <- rbindlist(signal)
+    upstream_signal <- signal[V2 %in% vp_idx][, c(1, 3, 4)]
+    downstream_signal <- signal[V1 %in% vp_idx][, 2:4]
   } else {
     # run for a region =========================================================
     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]
+    upstream_signal <- lapply(seq_along(explist), function(i) {
+      explist[[i]]$MAT[list(range_idx, vp_idx), 
+                       nomatch = 0][, list(V1, V3, exp = i)]
+    })
+    upstream_signal <- rbindlist(upstream_signal)
+    downstream_signal <- lapply(seq_along(explist), function(i) {
+      explist[[i]]$MAT[list(vp_idx, range_idx),
+                       nomatch = 0][, list(V2, V3, exp = i)]
+    })
+    downstream_signal <- rbindlist(downstream_signal)
   }
 
   signal <- rbind(upstream_signal, downstream_signal, use.names = FALSE)
-  colnames(signal) = c('V4','signal')
-  setkey(signal, 'V4')
+  colnames(signal) <- c('idx', 'signal', 'experiment')
+  signal <- signal[IDX, on = "idx==V4", nomatch = 0]
 
   # set basepairs --------------------------------------------------------------
-  bed_ranges <- exp$IDX[match(signal$V4, exp$IDX$V4)]
-  setkey(bed_ranges, 'V4')
-  signal <- signal[bed_ranges]
-  signal$mid = rowMeans(signal[,4:5])
+  signal$mid = rowMeans(signal[, 5:6])
 
   if( !is.null(xlim) ){
     signal <- signal[V1 == viewpoint[1,1]]
   }
 
 
   # output ---------------------------------------------------------------------
-  signal <- unique(signal[,c(3,6,2)])
-  colnames(signal) <- c('chromosome','mid','signal')
+  signal <- unique(signal[, c(4, 7, 2, 3)])
+  colnames(signal) <- c("chromosome", "mid", "signal", "experiment")
+  signal$experiment <- expnames[signal$experiment]
   signal <- list(data = signal)
 
   signal <- structure(signal,