fugi: Further Utilities for Genomic Interactions

#' Virtual 4C viewpoint
#'
#' Create a \linkS4class{GenomicInteractions} object with interactions originating at a given \dQuote{bait} viewpoint,
#' or a set of such viewpoints. 
#' This is similar to the idea of a virtual 4C experiment where you are interested in interactions with a specific region. 
#' 
#' @details
#' The object returned has the \dQuote{bait} as anchor one and the interacting regions as anchor two. 
#' By default this is genome wide.
#' If you only want to consider interactions within a certain distance around the bait, you can specify a region to consider.
#' 
#' Multiple baits can be given, e.g., to find all interactions around promoters.
#' 
#' You may want to visualise the resulting interactions in a genome browser - 
#' you can do this by creating coverage over anchor two of the object and exporting as a wig or bedgraph file.
#' 
#' @param x A \linkS4class{GenomicInteractions} object.
#' @param bait A \linkS4class{GRanges} object describing bait regions.
#' @param region An optional \linkS4class{GRanges} object specifying the region to look for bait interactions in.
#' @param ... Additional arguments to pass to \code{\link{findOverlaps}}.
#' 
#' @return A \linkS4class{GenomicInteractions} object.
#' 
#' @import GenomicRanges
#' @export
#' @examples
#' data(hic_example_data)
#' hic_example_data <- updateObject(hic_example_data)
#'
#' pos <- GRanges(seqnames="chr15", ranges=IRanges(start=59477709, end=59482708))
#' region <- GRanges(seqnames="chr15", ranges=IRanges(start=58980209, end=59980208))
#' vp <- viewPoint(hic_example_data, pos, region)
viewPoint = function(x, bait, region=NULL, ...) {
  hits <- list()  
  hits$one <- findOverlaps(x, bait, use.region = "first")
  hits$two <- findOverlaps(x, bait, use.region = "second")
  if (length(hits) == 0) return(GenomicInteractions())
  vp <- GenomicInteractions(anchor1 = bait[c(subjectHits(hits$one), 
                                          subjectHits(hits$two))],
                             anchor2 = c(anchorTwo(x)[queryHits(hits$one)], 
                                          anchorOne(x)[queryHits(hits$two)]),
                             counts = c(interactionCounts(x)[queryHits(hits$one)], 
                                      interactionCounts(x)[queryHits(hits$two)]))
  if (!is.null(region)) { 
    vp <- x[overlapsAny(anchorTwo(x), region, ...)] 
    }
  ord = order(vp)
  return(vp[ord])
}

#' Plot coverage around a virtual 4C viewpoint
#' 
#' Plots coverage of interactions around a given viewpoint. 
#' This function requires the output of `viewPoint()` as input. 
#' You should additionally specify the total region you wish to plot. 
#'
#' @param x A \linkS4class{GenomicInteractions} object generated by \code{\link{viewPoint}}.
#' @param region A \linkS4class{GRanges} object specifying the genomic region to plot.
#' @param ylab String containing the Y axis label.
#' @param xlab String containing the X axis label. 
#' By default this is the chromosome of the region that is being plotted.
#' @param ... Additional arguments to pass to \code{\link{plot}}.
#' 
#' @return A coverage track is plotted on the current graphics device.
#' A \linkS4class{RleList} object containing the result of \code{\link{coverage}} is invisibly returned.
#'
#' @examples
#' data(hic_example_data)
#' hic_example_data <- updateObject(hic_example_data)
#'
#' pos <- GRanges(seqnames="chr15", ranges=IRanges(start=59477709, end=59482708))
#' region <- GRanges(seqnames="chr15", ranges=IRanges(start=58980209, end=59980208))
#' vp <- viewPoint(hic_example_data, pos, region)
#' plotViewpoint(vp, region)
#'
#' @import GenomicRanges
#' @importFrom graphics plot
#' @export
#' 
plotViewpoint = function(x, region, ylab="Signal", xlab=NULL, ...) {
    if (length(region) > 1) stop("region must be a single range")
    x = x[overlapsAny(anchorTwo(x), region, type="within")]
    cov = as(coverage(anchorTwo(x), weight = interactionCounts(x))[region], "GRanges")
    points_x = c(start(cov), end(cov)) + start(region)
    points_y = rep.int(cov$score, 2)
    ord = order(points_x)
    
    if(is.null(xlab)){ 
    xlab <- as.character(seqnames(region))
    }
    
    plot(points_x[ord], points_y[ord], t="l", 
             ylab = ylab, xlab = xlab, ...)
    
    return(invisible(coverage(cov, weight = "score")))
}

#' Plot coverage around a set of virtual 4C viewpoints
#' 
#' Plots summarised coverage of interactions around a set of viewpoints, e.g., promoters. 
#' This function requires the output of \code{\link{viewPoint}} as input. 
#'
#' @param x A \linkS4class{GenomicInteractions} object generated by \code{\link{viewPoint}}.
#' @param left_dist Integer scalar specifying the distance to the left of each bait to consider, in bp.
#' @param right_dist Integer scalar specifying the distance to the right of each bait to consider, in bp.
#' @param fix One of \code{"center"}, \code{"start"}, \code{"end"}, to be passed to \code{\link{resize}}.
#' Interaction distances are calculated relative to this part of the bait. 
#' @param ylab String containing the Y axis label.
#' @param xlab String containing the X axis label. 
#' @param ... Additional arguments to \code{\link{plot}}.
#' 
#' @return A coverage track is plotted on the current graphics device.
#' A \linkS4class{RleList} object containing the aggregated \code{\link{coverage}} is invisibly returned.
#'
#' @examples
#' data(hic_example_data)
#' hic_example_data <- updateObject(hic_example_data)
#' 
#' pos <- GRanges(seqnames="chr15", ranges=IRanges(start=59477709, end=59482708))
#' region <- GRanges(seqnames="chr15", ranges=IRanges(start=58980209, end=59980208))
#' vp <- viewPoint(hic_example_data, pos, region)
#' plotAvgViewpoint(vp, left_dist = 1000000, right_dist = 100000)
#' @import GenomicRanges
#' @importFrom S4Vectors runValue
#' @export
plotAvgViewpoint = function(x, left_dist = 100000, right_dist = 100000, ylab="Average signal", 
                            xlab="Relative position", fix = "center",...) {
  
  cov <- .makeRelativeVP(x, fix = fix, left_dist = left_dist, right_dist = right_dist )
  
  points_x = c(start(cov), end(cov)) - left_dist
  points_y = rep.int(runValue(cov), 2)
  ord = order(points_x)
  
  p = plot(points_x[ord], points_y[ord], t="l", 
           ylab = ylab, xlab = xlab, ...)
  
  return(invisible(cov))
}

.makeRelativeVP <- function(GIObject, fix = "center", left_dist = 100000, right_dist = 100000){
  #make ranges relative to bait
  bait <- resize(anchorOne(GIObject), width = 1, fix = fix)
  ints <- ranges(anchorTwo(GIObject))
  ints <- GenomicRanges::shift(ints, shift = -(start(bait)))
  ints <- GenomicRanges::shift(ints, shift = left_dist)
  
  #make coverage and adjust to mean coverage per bait
  cov <- coverage(ints, weight = interactionCounts(GIObject), 
                  width = right_dist+left_dist)
  cov <- cov / length(unique(bait))
  
  return(cov)
  #get points to plot
  
}