R/ViewGuides.R

In RussBainer/gCrisprTools: Suite of Functions for Pooled Crispr Screen QC and Analysis

##' @title Assign Colors Based on the Position of a Value in a Distribution
##' @description This is a function to generate colors for plot elements on the basis of the position of a value within a distribution. 
##' Called internally by ct.viewGuides. 
##' @param exprs The value whose color is to be returned. 
##' @param rankedexprs A vector of values the length of cols that corresponds to the values in the distribution.
##' @param colors The vector of colors to be used as a reference. 
##' @return A value contained in cols. 
##' @keywords internal
##' @author Russell Bainer
ct.exprsColor <- function(exprs, rankedexprs, colors) {
    colors[which(abs(rankedexprs - exprs) == min(abs(rankedexprs - exprs)))]
}

##' @title Draw a horizontal line of a specified color.
##' @description This is a function called internally by ct.viewGuides to generate the color legend. End users should not use it. 
##' @param x,y Minimal coordinates to specify the line, which is drawn from the Y axis. 
##' @param color Guess!
##' @param width Line width. 
##' @return A line on an open device. 
##' @keywords internal
##' @author Russell Bainer
ct.drawFlat <- function(x, y, color, width = 1) {
    lines(c(0, x), c(y, y), col = color, lwd = width)
}

##' @title Draw a density color legend. 
##' @description This is a function called internally by ct.viewGuides to generate the color legend. End users should not use it. 
##' @param dens A density object. 
##' @param colorscale A vector of colors to draw behind the density. 
##' @return A color legend on the current graphics device. 
##' @keywords internal
##' @author Russell Bainer
ct.drawColorLegend <- function(dens, colorscale) {
    yrange <- (rep(max(abs(range(dens$x))), 2) * c(-1, 1))

    plot(NA, xlim = c(0, max(dens$y)), ylim = range(dens$x), xaxt = "n", xlab = "", ylab = "Log2 Normalized Reads", bty = "n")
    invisible(mapply(ct.drawFlat, x = rep(max(dens$y), length(colorscale)), y = seq(yrange[2], yrange[1], length.out = length(colorscale)), width = 2, color = colorscale))
    polygon(x = dens$y, y = dens$x, col = "black")
}

##' @title Generate a Plot of individual gRNA Pair Data in a Crispr Screen
##' @description This function generates a visualization of the effect estimates from a MArrayLM model result for all of the 
##' individual guides targeting a particular element, specified somewhere in the library annotation file. The estimated effect 
##' size and variance is plotted relative to zero for the specified contrast, with the color of the dot indicating the relative scale of the
##' of the guide intercept within the model framework, with warmer colors indicating lowly expressed guides. For comparison, the density of gRNA 
##' fold change estimates is privided in a pane on the right, with white lines indicating the exact levels of the individual guides. 
##' @param gene the name of the target element of interest, contained within the 'type' column of the annotation file.
##' @param fit An object of class MArrayLM containing, at minimum, an 'Amean' slot containing the guide level abundances, 
##' a 'coefficients' slot containing the effect estimates for each guide, and an 'stdev.unscaled' slot giving the coefficient standard Deviations. 
##' @param ann A data.frame object containing the gRNA annotations. 
##' At mimimum, it should have a column with the name specified by the \code{type} argument, containing the element targeted by each guide. 
##' @param type A character string indicating the column in ann containing the target of interest. 
##' @param contrast.term If a fit object with multiple coefficients is passed in, a string indiating the coefficient of interest.   
##' @param ylims An optional numeric vector of length 2 indicating the extremes of the y-axis scale. 
##' @return An image summarizing gRNA behavior within the specifed gene on the default device. 
##' @author Russell Bainer
##' @examples
##' data('fit')
##' data('ann')
##' ct.viewGuides('Target1633', fit, ann)
##' @export
ct.viewGuides <- function(gene, fit, ann, type = "geneSymbol", contrast.term = NULL, ylims = NULL) {
    current.graphic.params <- par(no.readonly = TRUE)
    on.exit(suppressWarnings(par(current.graphic.params)))

    par(xpd = FALSE)

    if (ncol(fit$coefficients) > 1) {
        if (is.null(contrast.term)) {
            stop("The fit object contains multiple coefficients. Please specify a contrast.term.")
        }
        fit <- ct.preprocessFit(fit, contrast.term)
    }



    # testing
    if (!methods::is(fit, "MArrayLM")) {
        stop(deparse(substitute(eset)), " is not an MArrayLM.")
    }

    # Find the gRNAs targeting the gene from the annotation, and order them
    options(warn = -1)
    ann <- ct.prepareAnnotation(ann, fit, controls = FALSE)
    options(warn = 0)

    if (!(sum(ann[, type] %in% gene))) {
        stop(gene, " is not present in the annotation file.")
    }

    grna.inx <- row.names(ann)[(ann[, type] %in% gene)]
    grna.inx <- sort(grna.inx)

    # Set up the color palette and the legend first:
    day3Density <- density(fit$coefficients)

    if (is.null(ylims)) {
        ylimit <- max(fit$coefficients)
        ylims <- range(c(min(0, (fit$coefficients[grna.inx, 1] - fit$stdev.unscaled[grna.inx, 1])), max(0, (fit$coefficients[grna.inx, 1] + fit$stdev.unscaled[grna.inx, 
            1]))))
    } else {
        if (any(length(na.omit(ylims)) != 2, !is.numeric(ylims), sum(is.infinite(ylims)) > 0)) {
            stop("ylims must be a numeric vector of length 2 containing finite values.")
        }
    }

    # Set up the guide colors based on their overall expression
    cols <- colorRampPalette(c("red", "orange", "white", "blue", "purple"))(1000)
    rankexprscolors <- seq(from = max(fit$Amean), to = min(fit$Amean), length.out = 1000)

    # Set up a couple of panes:
    layout(matrix(c(1, 2), ncol = 2, byrow = TRUE), widths = c(5, 0.5), heights = 3, respect = TRUE)
    par(mar = c(5, 4, 4, 2))

    plot(seq_len(length(grna.inx)), fit$coefficients[grna.inx, 1], ylab = "Log2(Fold Change)", main = gene, xlab = "", ylim = sort(ylims), col = vapply(fit$coefficients[grna.inx, 
        1], ct.exprsColor, character(1), rankedexprs = rankexprscolors, colors = cols), xaxt = "n", pch = 18)
    segments(seq_len(length(grna.inx)), fit$coefficients[grna.inx, 1], seq_len(length(grna.inx)), (fit$coefficients[grna.inx, 1] - fit$stdev.unscaled[grna.inx, 1]))
    segments(seq_len(length(grna.inx)), fit$coefficients[grna.inx, 1], seq_len(length(grna.inx)), (fit$coefficients[grna.inx, 1] + fit$stdev.unscaled[grna.inx, 1]))


    axis(1, at = seq_len(length(grna.inx)), labels = grna.inx, las = 2)
    abline(h = 0, lty = 4, col = "red")
    abline(v = (seq_len(length(grna.inx)) + 0.5), col = "grey")
    # title(main = )
    if (!is.null(contrast.term)) {
        mtext(contrast.term, side = 3, line = 0)
    }
    # Plop in the legend:
    par(mar = c(5, 1, 4, 1))

    ct.drawColorLegend(dens = day3Density, colorscale = cols)
    mapply(ct.drawFlat, x = rep(max(day3Density$y), length(grna.inx)), y = fit$coefficients[grna.inx, 1], width = 0.5, color = "black")
    lines(c(0, max(day3Density$y)), c(0, 0), col = "red", lty = "dashed")
    return(invisible(TRUE))
}