R/InteractionTrack-class.R

In LTLA/fugi: Further Utilities for Genomic Interactions

#' @title The InteractionTrack class
#'
#' @description
#' A class to hold chromatin interaction data for a specific genomic region.
#'
#' The InteractionTrack class is a specific \pkg{Gviz}-derived class for enabling the visualisation of chromatin interaction data. 
#' It allows interactions on a specified chromosome to be visualised by examining interactions between anchors as bezier curves. 
#' The object is instantiated and used in a similar fashion to standard \pkg{Gviz} tracks, 
#' e.g., plotted using the \code{\link{plotTracks}} function.
#' 
#' Several additional display parameters are defined for this class:
#' \describe{
#' \item{\code{plot.anchors}:}{Logical scalar indicating whether anchor elements are to be drawn.}
#' \item{\code{col.anchors.line}:}{String specifying the colour of border of these anchor elements.}
#' \item{\code{col.anchors.fill}:}{String specifying the colur to fill these anchor elements.}
#' \item{\code{plot.outside}:}{Logical scalar indicating whether interactions spanning outside of the window are to be plotted.}
#' \item{\code{col.outside}:}{String specifying the line colour to use for outside-spanning interactions.}
#' \item{\code{plot.trans}:}{Logical scalar indicating whether trans-chromosomal interactions are to be plotted.}
#' \item{\code{col.trans}:}{String specifying the line colour to use for trans-chromosomal interactions.}
#' \item{\code{col.interactions}:}{String specifying the colour of arcs representing interactions within the region of interest.}
#' \item{\code{col.interaction.types}:}{A named character vector of colours to use for each type of interaction,
#' where the name corresponds to the type of interaction (e.g., promoter-promoter, distal-promoter).}
#' \item{col.anchors.line.node.class}{A named character vector of colours to use for the border of each type of anchor element,
#' where the name corresponds to the anchor type (e.g., promoter).}
#' \item{col.anchors.fill.node.class}{A named character vector of colours to use for the fillof each type of anchor element,
#' where the name corresponds to the anchor type (e.g., promoter).}
#' }
#'
#' A number of options are available to convey quantitative information:
#' \itemize{
#' \item By default, the height of an arc representing an interaction is proportional to the number of supporting reads.
#' Instead of using the counts to define this, the height can be set to be proportional to either \code{fdr} or \code{p.value} using the \code{interaction.measure} display parameter. 
#' \item By changing the \code{interaction.dimension} to \code{"width"}, the line widths of each arc now represent the statistic supporting them, instead of the heights.
#' \item The heights of the arcs can be made to be proportional to the log10 of the supporting statistic by changing \code{interaction.dimension.transform} to \code{"log"}. 
#' }
#'
#' All of these parameters can be set with standard methods, i.e., \code{\link{displayPars<-}} on a InteractionTrack object.
#'
#' @import Gviz
#' @import grid
#' @importFrom stringr str_split
#'
#' @export InteractionTrack
setClass("InteractionTrack",
         contains=c("GdObject"),
         representation=representation(giobject="GenomicInteractions",
                                       variables="list",
                                       chromosome="character",
                                       stacking="character"),
         prototype=prototype(dp=DisplayPars(plot.anchors=TRUE,
                                            col.anchors.line = "lightblue",
                                            col.anchors.fill = "lightblue",
                                            col.anchors.line.node.class = c(),
                                            col.anchors.fill.node.class = c(),
                                            interaction.dimension="height",
                                            interaction.measure ="counts",
                                            interaction.dimension.transform = "default",
                                            col.interactions = "red",
                                            plot.outside = TRUE,
                                            col.outside = "red",
                                            plot.trans = FALSE,
                                            col.trans = "lightgray",
                                            col.interaction.types = c(),
                                            anchor.height = 0.05
                                            )))


setMethod("initialize", "InteractionTrack", function(.Object, giobject, chromosome, variables, ...) {
  .Object <- Gviz:::.updatePars(.Object, "InteractionTrack") #
  .Object@giobject <- giobject
  .Object@chromosome <- chromosome
  .Object@variables <- variables
  .Object <- callNextMethod(.Object, ...)
  return(.Object)
})

setMethod("start", "InteractionTrack", function(x){
  
  if(!is.null(x@variables$start)){
    tmp.start = min(c(start(anchorOne(x@giobject))[ seqnames(anchorOne(x@giobject)) ==  x@chromosome ],
        start(anchorTwo(x@giobject))[ seqnames(anchorTwo(x@giobject)) ==  x@chromosome ]))
  }else{
    return(x@variables$start)
  }
  return(tmp.start)
  } )

setMethod("end", "InteractionTrack", function(x){ 
  if(!is.null(x@variables$end)){
    tmp.end = max(c(end(anchorOne(x@giobject))[ seqnames(anchorOne(x@giobject)) ==  x@chromosome ],
                    end(anchorTwo(x@giobject))[ seqnames(anchorTwo(x@giobject)) ==  x@chromosome ]))
  }else{
    return(x@variables$start)
  }
  return(tmp.end)
})

setMethod("chromosome", "InteractionTrack", function(GdObject) GdObject@chromosome)

setMethod("subset", signature(x="InteractionTrack"), function(x, from, to, chromosome, ...){
  x@chromosome = chromosome
  x@giobject = subsetByFeatures(x@giobject, GRanges(chromosome, IRanges(from, to)))
  return(x)                                          
})

#' Constructor to create an InteractionTrack object
#'
#' Create an \linkS4class{InteractionTrack} object from an \linkS4class{GenomicInteractions} object to visualise a specified chromosome.
#'
#' @param x A \linkS4class{GenomicInteractions} object.
#' @param chromosome String specifying which chromosome to use to create a track.
#' @param name String specifying the name of the track.
#' @param start Integer scalar specifying the start location for the track.
#' @param end Integer scalar specifying the end location for the track.
#'
#' @return An \linkS4class{InteractionTrack} object.
#'
#' @examples 
#' library(Gviz)
#' anchor.one = GRanges(c("chr1", "chr1", "chr1", "chr1"), 
#'     IRanges(c(10, 20, 30, 20), width=5))
#' anchor.two = GRanges(c("chr1", "chr1", "chr1", "chr2"), 
#'     IRanges(c(100, 200, 300, 50), width=5))
#'
#' interaction_counts = sample(1:10, 4)
#' test <- GenomicInteractions(anchor.one, anchor.two, counts=interaction_counts)
#' metadata(test) <- list(description="this is a test", experiment_name="test")
#' interactions.track = InteractionTrack(name="Test", test, chromosome="chr1")                        
#' plotTracks(list(interactions.track), chromosome="chr1", from=0, to=500)
#' 
#' @export
InteractionTrack <- function(x, chromosome="", name=name(x), start=NULL, end=NULL){ 
  if(!(class(x)=="GenomicInteractions")){ stop("x must be a GenomicInteractions object")}

  all.seqnames <- seqlevels(.get_single_regions(x))
  if(chromosome !="" & !(chromosome %in% all.seqnames)){
    stop(paste("chromosome:", chromosome, "not found in seqlevels of the supplied GenomicInteractions object", sep=" "))
  }
  if(chromosome != ""){
    if(xor(is.null(start), is.null(end))){
      stop(paste("both start and end need to be specified"))
    }else if(is.null(start) & is.null(end)){
      x = x[as.vector(seqnames(anchorOne(x)) == chromosome | seqnames(anchorTwo(x)) == chromosome)]
    }else{
      x = subsetByFeatures(x, GRanges(chromosome, IRanges(start, end)))
    }
  }
	return(new("InteractionTrack", name=name, giobject=x, chromosome=chromosome, variables=list(chromosome=chromosome, start=start, end=end)))	
}

setMethod("drawGD", signature("InteractionTrack"), function(GdObject, minBase, maxBase, prepare=FALSE, subset=TRUE, ...){ 
  
  if(subset){
    GdObject <- subset(GdObject, chromosome=GdObject@chromosome, from=minBase, to=maxBase)
  }
  
  if(prepare){
      pushViewport(viewport(xscale=c(minBase, maxBase), yscale=c(0, 1) )) 
      popViewport(1)
      return(invisible(GdObject))
    }

  pushViewport(viewport(xscale=c(minBase, maxBase), yscale=c(0, 1)))
  if(length(GdObject@giobject)>0){
    plot.anchors = displayPars(GdObject, "plot.anchors")
    col.anchors.line = displayPars(GdObject, "col.anchors.line")
    col.anchors.fill = displayPars(GdObject, "col.anchors.fill")
    col.anchors.fill.node.class = displayPars(GdObject, "col.anchors.fill.node.class")
    col.anchors.line.node.class = displayPars(GdObject, "col.anchors.line.node.class")
    anchor.height = displayPars(GdObject, "anchor.height")
    interaction.dimension.transform = displayPars(GdObject, "interaction.dimension.transform")
  
    col.interactions = displayPars(GdObject, "col.interactions")
    plot.outside = displayPars(GdObject, "plot.outside")
    col.outside =  displayPars(GdObject, "col.outside")
    plot.trans = displayPars(GdObject, "plot.trans")
    col.trans = displayPars(GdObject, "col.trans")
    
    anchor_one_chr = as.character(seqnames(anchorOne(GdObject@giobject)))
    anchor_one_starts = start(anchorOne(GdObject@giobject))
    anchor_one_ends = end(anchorOne(GdObject@giobject))
    
    anchor_two_chr = as.character(seqnames(anchorTwo(GdObject@giobject)))
    anchor_two_starts = start(anchorTwo(GdObject@giobject))
    anchor_two_ends = end(anchorTwo(GdObject@giobject))
    
    anchor_one_midpoints = (anchor_one_starts + anchor_one_ends) / 2
    anchor_two_midpoints = (anchor_two_starts + anchor_two_ends) / 2
    
    if(displayPars(GdObject, "interaction.dimension")=="width" && displayPars(GdObject, "interaction.measure") == "counts"){
      lwds = 1+log(interactionCounts(GdObject@giobject))
    }else if(displayPars(GdObject, "interaction.dimension")=="width" && displayPars(GdObject, "interaction.measure") == "fdr"){
      lwds = 1-(GdObject@giobject$fdr)
    }else if(displayPars(GdObject, "interaction.dimension")=="width" && displayPars(GdObject, "interaction.measure") == "p.value"){
      lwds = 1-(GdObject@giobject$p.value)
    }else{
      lwds = rep(1, length(anchor_one_chr))
    }
    
    trans.indexes = which( anchor_one_chr != GdObject@chromosome | anchor_two_chr != GdObject@chromosome )
    outside.indexes = which( (anchor_one_chr == GdObject@chromosome & anchor_two_chr == GdObject@chromosome) & (anchor_one_midpoints <minBase | anchor_one_midpoints > maxBase | 
       anchor_two_midpoints < minBase | anchor_two_midpoints > maxBase) )
    inside.indexes = 1:length(anchor_one_chr)
    inside.indexes = which(!(inside.indexes %in% c(trans.indexes, outside.indexes)))
    
    cols = rep(col.interactions, length(anchor_one_chr))
    cols[trans.indexes] = col.trans
    cols[outside.indexes] = col.outside
    
    if(length(displayPars(GdObject, "col.interaction.types"))>0){
      col.interaction.types = displayPars(GdObject, "col.interaction.types")
      colour.map = str_split(names(col.interaction.types), "-")
      for(i in 1:length(colour.map)){
        cols[ isInteractionType(GdObject@giobject, colour.map[[i]][1], colour.map[[i]][2]) ] = col.interaction.types[i]
      }
    }
    
    if(plot.anchors){
      y.start = anchor.height
    }else{
      y.start = 0
    }
    
    if(plot.trans & length(trans.indexes) >0){
      for(i in trans.indexes){
        if( anchor_one_chr[i] != GdObject@chromosome ){
          if(anchor_two_midpoints[i] > ((minBase + maxBase)/2)){
            grid.curve(x1=anchor_two_midpoints[i], y1=y.start, x2=maxBase, y2=0.95,  
                       curvature = -1, default.units="native", gp=gpar(col=cols[i],lwd=lwds[i]))
          }else{
            grid.curve(x1=anchor_two_midpoints[i], y1=y.start, x2=minBase, y2=0.95,  
                       curvature = 1, default.units="native", gp=gpar(col=cols[i],lwd=lwds[i]))
          }      
        }else if( anchor_two_chr[i] != GdObject@chromosome){
          if(anchor_one_midpoints[i] > ((minBase + maxBase)/2)){
            grid.curve(x1=anchor_one_midpoints[i], y1=y.start, x2=maxBase, y2=0.95, 
                       curvature = -1, default.units="native", gp=gpar(col=cols[i],lwd=lwds[i]))
          }else{
            grid.curve(x1=anchor_one_midpoints[i], y1=y.start, x2=minBase, y2=0.95, 
                       curvature = 1, default.units="native", gp=gpar(col=cols[i],lwd=lwds[i]))
          }
        }
      }
    }
    
    if(displayPars(GdObject, "interaction.dimension")=="height" && displayPars(GdObject, "interaction.measure") == "counts"){
      if(interaction.dimension.transform == "log"){
        counts = 1 + log(interactionCounts(GdObject@giobject))
        curve.heights = anchor.height + (( 1 - anchor.height ) * (counts / max(counts)))
      }else{
        curve.heights = anchor.height + (( 1 - anchor.height ) * interactionCounts(GdObject@giobject)/max(interactionCounts(GdObject@giobject)))
      }
    }else if(displayPars(GdObject, "interaction.dimension")=="height" && displayPars(GdObject, "interaction.measure") == "fdr"){
      if(interaction.dimension.transform == "log"){
        fdr = GdObject@giobject$fdr
        fdr[is.infinite(log10(fdr))] = 2.2e-16
        fdr = -log10(fdr)
      }else{
        fdr = 1 - GdObject@giobject$fdr
      }
      curve.heights = anchor.height + (1-anchor.height) * (fdr / max(fdr))
    }else if(displayPars(GdObject, "interaction.dimension")=="height" && displayPars(GdObject, "interaction.measure") == "p.value"){
      if(interaction.dimension.transform == "log"){
        p.value = GdObject@giobject$p.value
        p.value[is.infinite(log10(p.value))] = 2.2e-16
        p.value = -log10(p.value)
      }else{
        p.value = 1 - GdObject@giobject$p.value
      }
      curve.heights = anchor.height + (1-anchor.height) * (p.value / max(p.value))
    }else{
      curve.heights = rep(1, length(anchor_one_chr))
    }
    
    if(plot.outside & length(outside.indexes)>0){
      xs = c()
      ys = c()
      for(i in outside.indexes){
        mdpt = (anchor_one_midpoints[i] + anchor_two_midpoints[i]) / 2
        xs = c(xs, c(anchor_one_midpoints[i], mdpt, anchor_two_midpoints[i]))
        ys = c(ys, c(y.start, curve.heights[i], y.start))
      } 
      grid.xspline(xs, ys, id.lengths=rep(3, length(outside.indexes)), 
                   shape=-1, default.units="native", 
                   gp=gpar(col=cols[outside.indexes], 
                           lwd=lwds[outside.indexes]))
    }
    if(length(inside.indexes)>0){
      xs = c()
      ys = c()
      for(i in inside.indexes){
        mdpt = (anchor_one_midpoints[i] + anchor_two_midpoints[i]) / 2
        xs = c(xs, c(anchor_one_midpoints[i], mdpt, anchor_two_midpoints[i]))
        ys = c(ys,  c(y.start, curve.heights[i], y.start))
      }
      grid.xspline(xs, ys, id.lengths=rep(3, length(inside.indexes)), 
                   shape=-1, default.units="native", 
                   gp=gpar(col=cols[inside.indexes], 
                           lwd=lwds[inside.indexes]))
    }  
    if(plot.anchors & length(GdObject@giobject) > 0){
      
      anchors = unique(c(anchorOne(GdObject@giobject), anchorTwo(GdObject@giobject)))
      
      col.anchors.fill = rep(col.anchors.fill, length(anchors))
      if(length(col.anchors.fill.node.class) > 0 & "node.class" %in% names(elementMetadata(anchors))){
        indexes = which(anchors$node.class %in% names(col.anchors.fill.node.class))
        col.anchors.fill[ indexes ] = col.anchors.fill.node.class[ anchors$node.class[indexes] ]
      } 
        
      col.anchors.line = rep(col.anchors.line, length(anchors))
      if(length(col.anchors.line.node.class) > 0 & "node.class" %in% names(elementMetadata(anchors))){
        indexes = which(anchors$node.class %in% names(col.anchors.line.node.class))
        col.anchors.line[ indexes ] = col.anchors.line.node.class[ anchors$node.class[indexes] ]
      } 
      
      xs = c()
      ys = c()
      for(i in 1:length(anchors)){
        xs = c(xs, c(start(anchors[i]), end(anchors[i]), end(anchors[i]), start(anchors[i])))
        ys = c(ys, c(0, 0, anchor.height, anchor.height))
      }
      grid.polygon(xs, ys, id.lengths=rep(4, length(anchors)), 
                   default.units="native", 
                   gp=gpar(col=col.anchors.line, 
                           fill= col.anchors.fill))
    }
  }
  popViewport(1)
  return(invisible(GdObject))
})

#' Available display parameters 
#'
#' Query the default display parameters for an \linkS4class{InteractionTrack} object.
#' 
#' @param class A valid track object class name, or the object itself, in which case the class name is determined from it.
#'
#' @details
#' This function provides the same functionality as \code{Gviz::availableDisplayPars} and allows the user to display the default display parameters for the \code{InteractionTrack} class. 
#' If the class of the track is not \code{"InteractionTrack"} then the function calls the \code{\link{availableDisplayPars}} method in \pkg{Gviz}.
#' 
#' @return A named list of the default display parameters.
#' 
#' @examples
#' availableDisplayPars("InteractionTrack")
#' 
#' @importFrom Gviz availableDisplayPars
#' @export
availableDisplayPars = function(class){
  if(!is.character(class))
    class <- class(class)
  
  if(class=="InteractionTrack"){
    parents <- names(getClassDef(class)@contains)
    
    pars =try(sapply(c(parents, "InteractionTrack"), function(x) 
      as(getClassDef(x)@prototype@dp, "list"), simplify=FALSE), silent=TRUE)
    finalPars <- list()
    inherited <- list()
    for(p in names(pars))
    {
      finalPars[names(pars[[p]])] <- pars[[p]]
      inherited[names(pars[[p]])] <- p
    }
    finalPars <- finalPars[order(names(finalPars))]
    inherited <- inherited[order(names(inherited))]
    return(new("InferredDisplayPars", name=class, 
               inheritance=unlist(inherited), finalPars)) 
  }else{
    Gviz::availableDisplayPars(class)
  }
}