R/plotBedgraph.R

In dphansti/Sushi: Tools for visualizing genomics data

#' plots data stored in bed file format
#' 
#'
#' @param signal signal track data to be plotted (in bedgraph format)
#' @param chrom chromosome of region to be plotted
#' @param chromstart start position
#' @param chromend end position
#' @param range y-range to plpt ( c(min,max) )
#' @param color color of signal track  
#' @param lwd color of line outlining signal track.  (only valid if linecol is not NA)
#' @param linecolor color of line outlining signal track.  use NA for no outline
#' @param addscale TRUE/FALSE whether to add a y-axis
#' @param overlay TRUE / FALSE whether this data should be plotted on top of an existing plot
#' @param rescaleoverlay TRUE/FALSE whether the new plot shold be rescaled based on the maximum value to match the existing plot (only valid when overlay is set to 'TRUE')   
#' @param transparency Value between 0 and 1 indication the degree of transparency of the plot
#' @param flip TRUE/FALSE whether the plot should be flipped over the x-axis
#' @param xaxt A character which specifies the x axis type.  See \code{\link{par}}
#' @param yaxt A character which specifies the y axis type.  See \code{\link{par}}
#' @param xlab Label for the x-axis
#' @param ylab Label for the y-axis
#' @param xaxs Must be set to 'i' for appropriate integration into Sushi plots.  See \code{\link{par}}
#' @param yaxs Must be set to 'i' for appropriate integration into Sushi plots.  See \code{\link{par}}
#' plottype
#' @param bty A character string which determined the type of box which is drawn about plots.  See \code{\link{par}}
#' @param ymax fraction of max y value to set as height of plot.
#' @param colorbycol palette to use to shade the signal track plot.  Only applicable when overlay is set to FALSE.
#' @param ... values to be passed to \code{\link{plot}}
#' @export
#' @examples
#'
#' data(Sushi_ChIPSeq_CTCF.bedgraph)
#' data(Sushi_DNaseI.bedgraph)
#'
#' chrom            = "chr11"
#' chromstart       = 1955000
#' chromend         = 1965000
#' 
#' plotBedgraph(Sushi_ChIPSeq_CTCF.bedgraph,chrom,chromstart,chromend,transparency=.50,flip=FALSE,color="blue",linecol="blue")
#' plotBedgraph(Sushi_DNaseI.bedgraph,chrom,chromstart,chromend,transparency=.50,flip=FALSE,color="#E5001B",linecol="#E5001B",overlay=TRUE,rescaleoverlay=TRUE)
#' labelgenome(chrom,chromstart,chromend,side=1,scipen=20,n=3,line=.18,chromline=.5,scaleline=0.5,scale="Mb")
#' 
#' transparency = 0.5
#' col1 = col2rgb("blue")
#' finalcolor1 = rgb(col1[1],col1[2],col1[3],alpha=transparency * 255,maxColorValue = 255)
#' col2 = col2rgb("#E5001B")
#' finalcolor2 = rgb(col2[1],col2[2],col2[3],alpha=transparency * 255,maxColorValue = 255)
#' 
#' legend("topright",inset=0.025,legend=c("DnaseI","ChIP-seq (CTCF)"),fill=c(finalcolor1,finalcolor2),border=c("blue","#E5001B"),text.font=2,cex=0.75)
plotBedgraph <-
  function(signal,chrom,chromstart,chromend,range=NULL,color=SushiColors(2)(2)[1],
           lwd=1,linecolor=NA,addscale=FALSE,overlay=FALSE,rescaleoverlay=FALSE,transparency=1.0,
           flip=FALSE, xaxt='none',yaxt='none',xlab="",ylab="",xaxs="i",yaxs="i",bty='n',ymax=1.04,
           colorbycol=NULL,...)
  {
    if (overlay == TRUE)
    {
      colorbycol = NULL
    }
    
    if(is.na(linecolor ) == TRUE)
    {
      linecolor = color
    }
    
    # ensure that the chromosome is a character
    signal[,1] = as.character(signal[,1])
    
    # filter for desired region
    signaltrack = signal[which(signal[,1] == chrom & ((signal[,2] > chromstart & signal[,2] < chromend) |  (signal[,3] > chromstart & signal[,3] < chromend))),(2:4)]
    
    # remove any duplicate rows
    signaltrack = signaltrack[!duplicated(signaltrack),]
    
    # exit if overlay is TRUE and there isn't enough data
    if (overlay ==TRUE && nrow(signaltrack) < 2)
    {
      return ("not enough data within range to plot")
    }
    
    # exit if overlay is FALSE and there isn't enough data
    if (nrow(signaltrack) < 2)
    {
      if (is.null(range) == TRUE)
      {
        range = c(0,1)
      }
      
      # make blank plot
      plot(0,0,xlim=c(chromstart,chromend),type='n',ylim=range,xaxt=xaxt,yaxt=yaxt,ylab=ylab,xaxs=xaxs,yaxs=yaxs,bty=bty,xlab=xlab,...) 
      return ("not enough data within range to plot")
    }
    
    # downsample for plotting
    while (nrow(signaltrack) > 8000)
    {
      # downsample for plotting if neccesary
      if (nrow(signaltrack) %% 2 != 0)
      {
        signaltrack = signaltrack[1:(nrow(signaltrack)-1),]
      }
      
      starts = signaltrack[seq(1, nrow(signaltrack), 2),1]
      stops  = signaltrack[seq(2, nrow(signaltrack), 2),2]
      meanval = apply(cbind(signaltrack[seq(1, nrow(signaltrack), 2),3],signaltrack[seq(2, nrow(signaltrack), 2),3]), 1, mean)
      signaltrack = cbind(starts,stops,meanval)
    }
    
    # add col names
    names(signaltrack)[(1:3)]    = c("V1","V2","V3")
    
    # make linking regions if neccesary
    linkingregions = cbind(signaltrack[1:(nrow(signaltrack)-1),2], signaltrack[2:nrow(signaltrack),1])
    linkingregions = matrix(linkingregions[which(linkingregions[,1] != linkingregions[,2]),],ncol=2)
    
    if (nrow(linkingregions) > 0)
    {
      linkingregions = cbind(linkingregions,0)
      
      # make col names the same
      names(linkingregions)[(1:3)] = c("V1","V2","V3")
      
      # add linking regions to signaltrack
      signaltrack = rbind(signaltrack,linkingregions)
    }
    
    # sort data
    signaltrack = signaltrack[order(signaltrack[,1]),]
    
    # convert two columns to one
    signaltrack = cbind(as.vector(t(signaltrack[,c(1,2)])),as.vector(t(signaltrack[,c(3,3)])))
    
    # add slighltly negative vaue to both ends to ensure proper polygon plotting
    signaltrack = rbind(c(min(signaltrack[,1]),-.00001),signaltrack)
    signaltrack = rbind(signaltrack, c(max(signaltrack[,1]),-.00001))
    
    if (flip == TRUE)
    {
      signaltrack[,2] = signaltrack[,2]*-1
    }
    
    # determine the y-limits
    if (is.null(range) == TRUE)
    {
      range = c(0,ymax*max(signaltrack[,2]))
      if (flip == TRUE)
      {
        range = c(ymax*min(signaltrack[,2]),0)
      }
    }
    
    if (overlay == FALSE)
    {
      # make blank plot
      plot(signaltrack,xlim=c(chromstart,chromend),type='n',ylim=range,xaxt=xaxt,yaxt=yaxt,ylab=ylab,xaxs=xaxs,yaxs=yaxs,bty=bty,xlab=xlab,...) 
    }
    
    # rescale the overlay plot for comparative purposes
    if (rescaleoverlay == TRUE)
    {
      if (flip == FALSE)
      {
        signaltrack[,2] =   par('usr')[4] * signaltrack[,2] / max(abs(signaltrack[,2]) )
      }
      if (flip == TRUE)
      {
        
        signaltrack[,2] =   abs(par('usr')[3]) * signaltrack[,2] / max(abs(signaltrack[,2]) )
      }
    }
    
    # set the transparency
    rgbcol = col2rgb(color)
    finalcolor = rgb(rgbcol[1],rgbcol[2],rgbcol[3],alpha=transparency * 255,maxColorValue = 255)
    
    if (is.null(colorbycol) == FALSE)
    {
      # add the gradient to the background
      if (is.null(colorbycol) == FALSE)
      {
        plotlef = par('usr')[1]
        plotrig = par('usr')[2]
        plotbot = par('usr')[3]
        plottop = par('usr')[4]
        
        bgcol = maptocolors((1:100),colorbycol)
        
        if (flip == FALSE)
        {
          tops = seq(plotbot,plottop,length.out=101)[2:101]  
          bots =  seq(plotbot,plottop,length.out=101)[1:100]
        }
        if (flip == TRUE)
        {
          bots =  rev(seq(plotbot,plottop,length.out=101)[2:101])
          tops =  rev(seq(plotbot,plottop,length.out=101)[1:100] )
        }
        
        for (i in (3:(nrow(signaltrack)-1)))
        {
          if (flip == FALSE)
          {
            ybots = bots[which(bots <= signaltrack[i,2])]
            ytops = tops[which(tops <  signaltrack[i,2])]
          }
          if (flip == TRUE)
          {
            ybots = bots[which(bots >= signaltrack[i,2])]
            ytops = tops[which(tops >  signaltrack[i,2])]
          }
          ytops = c(ytops,signaltrack[i,2] )
          xleft = rep(signaltrack[i-1,1],length(ytops))
          xrigh = rep(signaltrack[i,1],length(ytops))
          xybgcol = bgcol[1:length(ytops)]
          
          rect(xleft=xleft, ybottom=ybots, xright=xrigh, ytop=ytops,col= xybgcol,border=bgcol,lwd=lwd) 
        } 
      }
    }
    
    if (is.null(colorbycol) == TRUE)
    {
      # plot the signal track
      polygon(signaltrack,col=finalcolor,border=linecolor,lwd=lwd)
    }
    
    # add scale to upper right corner
    if (addscale == TRUE)
    {
      mtext(paste(range[1],range[2],sep="-"),side=3,font=1,adj=1.00,line=-1)
    }
  }