R/plotRDI.R

In rdi: Repertoire Dissimilarity Index

#' RDI Axis annotation function
#' 
#' @name rdiAxis
#' @description This function takes a RDI model, as generated by \link{rdiModel}, 
#'   and adds an axis with annotations in the fold change space.
#' 
#' @param model     A model object, as generated by \link{rdiModel}.
#' @param side      The side the axis will be added to. (1 - bottom; 2 - left; 3 - top;
#'                  4 - right). Default is 2.
#' @param at        The points at which the tick marks are drawn. By default, 
#'                  tickmarks are placed at 'round' fold/percent change values
#'                  using the "\link[base]{pretty}" breakpoints function. This may
#'                  not be ideal if log-RDI values are being plotted.
#' @param ...       Additional parameters to pass to \link[graphics]{axis}
#'
#' @details This function is designed to replace the default axes generated by a plot 
#' function. Instead of annotating the true RDI value, \code{rdiAxis} will estimate
#' the "true difference" values at various points within the plotting region, and will
#' annotate the axis with those estimates. 
#' 
#' It is worth noting that although the RDI value can range below \code{rdiModel}'s 
#' estimate for "identical" repertoires, no negative true difference values will be 
#' annotated, as these values do not make sense. 
#' 
#' @seealso \link{rdiModel}, \link{rdiLadder}, \link{plotRDIladder}
#'
#' @examples 
#' 
#' #create genes
#' genes = sample(letters, 10000, replace=TRUE)
#' #create sequence annotations
#' seqAnnot = data.frame(donor = sample(1:10, 10000, replace=TRUE))
#' #calculate RDI
#' d = rdi(genes, seqAnnot)
#' 
#' ##create a "baseVect" with the same probability as our features
#' ##since we sampled uniformly, the base vector has equal probability
#' baseVect = rep(1/length(letters),length(letters))
#' 
#' ##generate an RDI model
#' m = rdiModel(attr(d, "nseq"), baseVects=baseVect)
#' 
#' ##convert RDI to lfc
#' td = convertRDI(d,models=m)$pred
#' 
#' par(mar=c(4,4,1,4),las=1,mgp=c(3,0.5,0))
#' plot(td,d, ylab="RDI", xlab="LFC")
#' 
#' ##now add "true difference" axis annotation to the right side of the plot
#' rdiAxis(m, side=4)
#'
#' @export
rdiAxis <- function(model, side=2, at=NULL, ...){
  
  ##if the axis labels aren't specified, calculate them.
  if(is.null(at)){
    ##Figure out the axis range. 
    usr = par("usr")
    min = usr[ifelse(side %in% c(1,3), 1,3)]
    max = usr[ifelse(side %in% c(1,3), 2,4)]
    
    ##convert range to units
    fc.range = predict(model$rev.fit,c(min,max))$y
    fc.range[1] = max(0,fc.range[1])   ##the min FC can't go below 0
    at = pretty(fc.range)
  }
  
  ##calculate the true locations of "at"
  at.true = predict(model$fit, at)$y
  
  axis(side, at=at.true, labels=at,...)
}

#' RDI ladder plotting function
#' 
#' @name plotRDIladder
#' @description function for adding a pre-computed RDI ladder onto a plot
#' 
#' @param ladder        the ladder object to add, as created by \link{rdiLadder}
#' @param side          integer; value between 1 and 4 indicating where the ladder will be 
#'                      added. 1 - bottom, 2 - left, 3 - top, 4 - right.
#' @param toPlot        logical vector; which ladders should be plotted? By default, 
#'                      ladders that are significantly overlapped by their neighbor and
#'                      those that are majority outside the plotting region are removed.
#' @param labelLadder   logical; if \code{TRUE}, each curve will be annotated with the ladder name  
#' @param add           logical; if \code{TRUE}, the ladder will be added to the current plot
#' @param cex           character expansion for ladder labels.               
#' @param lineCol      the colors to be used for the ladder border. If the length of \code{col} 
#'                      is less than the length of ladder, \code{col} will be recycled.
#' @param fillCol      the colors to be used to fill the ladder. If the length of \code{col} 
#'                      is less than the length of ladder, \code{col} will be recycled.
#' 
#' @details 
#' This function is used in conjunction with \link{rdiLadder} to add a useful 
#' annotation to any plot containing RDI values. 
#' 
#' Because RDI values vary according to the number of genes and size of the repertoires,
#' they are not useful as numbers by themselves. Instead, it is useful to compare them 
#' with estimates of the true difference between the two repertoires. This function
#' adds a series of density curves along one side of a standard plotting region, each
#' one representing the most likely RDI values between two repertoires that vary by a set
#' amount. 
#' 
#' By default, not all density curves from the \code{ladder} parameter are plotted. 
#' Instead, the function intelligently chooses which ladders to plot based on the amount
#' of overlap between neighboring ladders. If a ladder is significantly overlapped by the
#' ladder below it, then the ladder will not be plotted. In addition, if the mean of a 
#' ladder is outside the main plotting region, it will be dropped. In order to control this 
#' behavior, you can directly specify which ladders are plotted using the \code{toPlot} 
#' parameter.
#' 
#' @return 
#' Invisibly returns the location of the ladder (if side 1 or 3, the y location; 
#' otherwise, the x location).
#' 
#' @seealso  \link{rdiLadder}, \link{rdiModel}, \link{rdiAxis}
#' 
#' @export
plotRDIladder = function(ladder, side=4, toPlot=NULL,
                         labelLadder=TRUE, add=TRUE, cex=0.7, 
                         lineCol=NULL,fillCol="#AAAAAA"){
  ##if lineCol not specified, use par("col")
  if(is.null(lineCol)){lineCol = par("col")}
  
  ##make cols the same length as ladder
  if(length(lineCol)<length(ladder)){lineCol=rep(lineCol, length.out=length(ladder))}
  if(length(fillCol)<length(ladder)){fillCol=rep(fillCol, length.out=length(ladder))}
  
  ##plotting params
  srt = c(0,-90,0,90)[side]
  sign = ifelse(side %in% c(1,2), 1,-1)
  
  ##figure out the sides
  horiz = FALSE
  if(side %in% c(1,3)){horiz = TRUE; side = side+1}
  
  if(!add){
    ladder.lims = unlist(lapply(ladder, function(m){
      qnorm(c(0.001,0.999), mean=m["mean"], sd=m["sd"])
    }))
    ladder.rng = range(ladder.lims)
    xlim=c(0,1);ylim=ladder.rng
    if(horiz){xlim=ladder.rng;ylim=c(0,1);}
    plot(0, type="n", yaxt=ifelse(horiz,"n","s"),xaxt=ifelse(horiz,"s","n"),
         xlim=xlim, ylim=ylim, xlab=ifelse(horiz,"RDI",""), ylab=ifelse(horiz,"","RDI")
        )
  }
  
  
  ##figure out the label size
  if(horiz){
    lws = strwidth(names(ladder), cex=cex)
    lh = strheight(names(ladder), cex=cex)
  }else{
    lws = strwidth.rot(names(ladder),cex=cex)
    lh = strheight.rot(names(ladder)[1],cex=cex)
  }
  
  ##decide which ladders to plot automatically
  if(!is.null(toPlot)){
    if(length(toPlot)!=length(ladder)){stop("length of toPlot must match ladder.")}
    plotLadder=toPlot
  }else{
    plotLadder = rep(T, length(ladder))
    if(length(ladder)>1){
      for(i in 2:length(ladder)){
        mu = ladder[[i]]["mean"]; s = ladder[[i]]["sd"]
        mu_p = ladder[[i-1]]["mean"]; s_p = ladder[[i-1]]["sd"]
        ##set of conditions for plotting
        plotLadder[i] =
          ##check if it overlaps the previous one too much
          !(plotLadder[i-1] && (mu-lws[i]/2)<(mu_p+s_p)) &&
          ##check if the top half is outside the plot range
          (mu+lws[i])<par("usr")[ifelse(horiz, 2,4)]
      }
    }
  }
  
  ##figure out the location of the edge we're working with
  yloc.ladder = par("usr")[2*horiz + 1 + (side==4)]
  ##put the ladder just inside that edge
  yloc.ladder = yloc.ladder + sign *  ifelse(horiz, strwidth.rot("|"), strwidth("|"))
  
  ##counting down so that the lower one is always on top
  for(i in rev(which(plotLadder))){  
    mu = ladder[[i]]["mean"]; s = ladder[[i]]["sd"]
    
    ##calculate the normal curve
    rng = qnorm(c(0.001,0.999), mean=mu, sd=s)
    
    x = seq(rng[1],rng[2],length.out=30)
    d = dnorm(x, mean=mu, sd=s)
    d = d/max(d) * lh * 1.5  ##set the height of the curve to 1.5x the label height
    ##figure out y location from heights
    d = yloc.ladder + sign * d 
    
    ##if vertical, switch x and y
    if(!horiz){tmp = d; d = x; x = tmp;}
    
    ##plot the curve
    polygon( x,d, col = fillCol[i], lty = 0  )
    lines(   x,d, col = lineCol[i], lwd = 1.2)
    #lines(range(d),c(mu,mu))
    
    ##label the ladder
    if(labelLadder){
      x = ifelse(horiz, mu, yloc.ladder); y=ifelse(horiz, yloc.ladder,mu)
      adj = c(0.5,-0.1); if(horiz && side==4){adj[2]=1.1}
      text(x,y,labels = names(ladder)[i],adj=adj,cex=cex,srt=srt)
    }
  }
  invisible(yloc.ladder)
}