R/Plotting.R

In fossbert/binilib: Carlo's little helper

#' Quick changes to the plotting environment
#'
#' @param a number of rows (passed to mfrow)
#' @param b number of columns (passed to mfrow)
#' @param brewer_n number of colours to retrieve from RColorBrewer palette
#' @param brewer_name RColorBrewer palette
#' @param ... further arguments passed to par
#' @return By default only a new colour palette
#' @export
mypar <- function(a = 1, b = 1, brewer_n = 9, brewer_name = "Set1", ...) {

    graphics::par(mfrow = c(a, b), ...)
    grDevices::palette(RColorBrewer::brewer.pal(brewer_n, brewer_name))
    tmp <- ggplot2::theme_bw() +
        ggplot2::theme(text = element_text(size = 8, family = 'Arial', color = 'black'),
                       plot.title = element_text(size = 9),
                       axis.text = element_text(size = 9))
    ggplot2::theme_set(tmp)
}

#' Plot quantiles per column for a numeric matrix
#'
#' @param eset numeric matrix typically containing gene expression data
#' @param quants integer of number of top MR to consider
#' @param brewer_name RColorBrewer palette, defaults to Dark2
#' @param ... further arguments passed to matplot
#' @return graph with the respective quantiles across the sample space
#' @export

kaboxplot <- function(eset,
                      quants = c(0.05, 0.25, 0.5, 0.75, 0.95),
                      brewer_name = 'Dark2',
                      ylab = 'Normalized expression',
                      xlab = 'Samples', ...) {

    nq <- length(quants)
    qs <- t(apply(eset, 2, quantile, probs = quants, na.rm = TRUE))
    info <- colMeans(qs)

    oldmar <- par()$mar
    par(mar = rep(4.1, 4), las = 1)
    matplot(qs,
            type = "l",
            lty = 1,
            col = RColorBrewer::brewer.pal(nq, brewer_name),
            ylab = ylab,
            xlab = xlab,
            ...)
    axis(side = 4, at = info, labels = names(info)) # info on quantiles
    par(mar = oldmar, las = 0)
}


#' Add a custom color gradient to a base R plot
#'
#' @param col a chracter vector as returned from colorRampPalette
#' @param level a factor generated by cutting a continuous variable into many pieces (e.g. 100)
#' @param side where to draw the color gradient, default is 4, i.e. right
#' @return added a color gradient to existing plot
#' @export
legend_col <- function(col, level, side = 4){

    opar <- par

    n <- length(col)

    bx <- par("usr")

    box.cx <- c(bx[2] + (bx[2] - bx[1]) / 1000,
                bx[2] + (bx[2] - bx[1]) / 1000 + (bx[2] - bx[1]) / 50)
    box.cy <- c(bx[3], bx[3])
    box.sy <- (bx[4] - bx[3]) / n

    xx <- rep(box.cx, each = 2)

    par(xpd = TRUE)
    for(i in 1:n){

        yy <- c(box.cy[1] + (box.sy * (i - 1)),
                box.cy[1] + (box.sy * (i)),
                box.cy[1] + (box.sy * (i)),
                box.cy[1] + (box.sy * (i - 1)))
        polygon(xx, yy, col = col[i], border = col[i])

    }
    par(new = TRUE)
    plot(0, 0, type = "n",
         ylim = c(min(level), max(level)),
         yaxt = "n", ylab = "",
         xaxt = "n", xlab = "",
         frame.plot = FALSE)
    axis(side = side, las = 2, tick = FALSE, line = .25)
    par <- opar
}


#' Split violin plot
#'
#' This type of plot takes a numeric vector and splits it twice. Once using a first categorical variable with
#' any number of levels and then and then a second time using a second categorical variable with
#' exactly TWO levels. A density curve will be estimated and side-by-side violin plots will be produced.
#'
#' @param x a numeric variable to be split in two ways
#' @param s1 categorical variable used for the first split
#' @param s2 categorical variable used for the second split
#' @param cols vector of lenght 2 indicating colors in R (hex, integer or plane name)
#' @param ylb default for y-axis, could be changed if you know the nature of x (e.g. Z-score)
#' @param rug logical, whether to plot the data points as a rug
#' @param legpos character, indicating the position of the legend
#' @param rotate_xlabels integer, if not NULL, will be used as degree rotation of x labels
#' @return a plot will be send to the graphics output
#' @export

split_violin <- function(x, s1, s2,
                         cols = c('cornflowerblue', 'salmon'),
                         ylb = '',
                         rug = FALSE,
                         legpos = 'topleft',
                         rotate_xlabels = NULL,
                         ...) {

    omar <- par()$mar
    omgp <- par()$mgp

    s1 <- factor(s1)
    s2 <- factor(s2)

    # check whether input is legit
    if (length(x) != length(s1) || length(x) != length(s2)) {
        stop("Careful, length of input vectors don't match up")
    }
    if (length(levels(s2)) != 2) stop("Sorry friend, s2 MUST have 2 levels")

    # split both x and s2 by s1
    l1 <- split(x = x, f = s1)
    splitlist <- split(s2, f = s1)

    # second split and density estimates
    l2 <- lapply(seq_along(l1), function(i) {

        tmp <- split(l1[[i]], splitlist[[i]])

        lapply(tmp, function(j) {
            tmp2 <- density(j, na.rm = TRUE)
            list(x = tmp2$x, y = tmp2$y,
                 mode_x = tmp2$x[which.max(tmp2$y)],
                 mode_y = tmp2$y[which.max(tmp2$y)],
                 original_x = j)
        })

    })
    names(l2) <- names(l1)

    # prepare plotting - careful: x is y and y is x

    # find max sum of consecutive modes in plot
    ml <- vapply(l2, function(i) i[[1]][['mode_y']], FUN.VALUE = numeric(1))
    mr <- vapply(l2, function(i) i[[2]][['mode_y']], FUN.VALUE = numeric(1))
    xm <- max(ml[-1] + mr[-length(mr)], na.rm = TRUE) + .25

    at <- seq(from = 1, by = xm, length.out = length(l2))
    xr <- range(density(x)$x)
    xr_simple <- seq(round(xr[1]), round(xr[2]), by = 1)

    if (ml[1] > 0.5 || mr[length(mr)] > 0.5) {
        xside <- max(c(ml[1], mr[length(mr)]), na.rm = TRUE)
        yr <- c(min(at)-xside, max(at)+xside)
    } else {
        yr <- c(min(at)-0.5, max(at)+0.5)
    }

    # define plotting environment, space below x depends on length of x labels
    par(mar = c(2.1, 2.1, 0.1, 0.1))
    plot(0,
         type='n',
         axes=FALSE,
         ylab="",
         xlab="",
         xlim = yr,
         ylim = xr,
         ...)

    # side-by-side violin plots
    for (i in seq_along(l2)){

        tmp <- l2[[i]]

        # left hand side
        x1 <- tmp[[1]][['x']]
        y1 <- tmp[[1]][['y']]
        my1 <- tmp[[1]][['mode_y']]
        mx1 <- tmp[[1]][['mode_x']]

        polygon(x = c(at[i]-y1, rep(at[i], length(y1))),
                y = c(x1, rev(x1)), col = cols[1])
        lines(x = c(at[i]-my1, at[i]), y = rep(mx1, 2))

        # right hand side
        x2 <- tmp[[2]][['x']]
        y2 <- tmp[[2]][['y']]
        my2 <- tmp[[2]][['mode_y']]
        mx2 <- tmp[[2]][['mode_x']]
        polygon(x = c(at[i]+y2, rep(at[i], length(y2))),
                y = c(x2, rev(x2)), col = cols[2])
        lines(x = c(at[i]+my2, at[i]), y = rep(mx2, 2))

        if (rug) {
           # yrug <- quantile(c(y1, y2), .25)
            yrug <- .05
            orgs1 <- tmp[[1]][['original_x']]
            for (k in seq_along(orgs1)) {
                lines(x = c(at[i]-yrug, at[i]), y = rep(orgs1[k], 2))
            }
            orgs2 <- tmp[[2]][['original_x']]
            for (k in seq_along(orgs2)) {
                lines(x = c(at[i]+yrug, at[i]), y = rep(orgs2[k], 2))
            }
        }
    }

    if (!is.null(rotate_xlabels)){
        text(x = at, y = par("usr")[3],
             labels = names(l2),
             srt = rotate_xlabels,
             adj = c(1,1),
             xpd = TRUE,
             ...)
    } else {
        mtext(text = names(l2), side = 1, at = at, line = 0, ...)
    }

    axis(side = 2, tck = -0.01, mgp = c(1,.3,0), las = 1, ...)
    mtext(text = ylb, side=2, line=1, ...)
    box()
    legend(legpos, legend = levels(s2), pch = 15, col = cols, ...)

    par(mar = omar, mgp = omgp, las = 0) # restore old settings

}


# Plotting functions
plothm. <- function(x, color=c("royalblue","firebrick2"), gama=1, grid=T, scmax=0, box=TRUE, ...) {
    coli <- colorScale(x=filterRowMatrix(x, nrow(x):1), color=color, gama=gama, scmax=scmax)
    image(1:ncol(x), 1:nrow(x), t(matrix(1:(ncol(x)*nrow(x)), nrow(x), ncol(x))), col=coli, ylab="", xlab="", axes=F, ...)
    if (box) box()
    if (grid) grid(ncol(x), nrow(x), col="black", lty=1)
}

#' colorScale
#'
#' This function generates a color scale
#'
#' @param x Vector or matrix of numeric values
#' @param color Vector of character strings indicating the colors for the scale. Up to three colors can be defined. While is used for the missing color
#' @param gama Number indicating the gama transformation
#' @param alpha Number between 0 and 1 indicating the transparency of the color (1 for absolute color)
#' @param scmax Number indicating the maximum value for the scale
#' @param nacol Character string indicating the color for missing values
#' @return Vector of colors
#' @export
colorScale <- function(x, color=c("royalblue","firebrick2"), gama=1, alpha=1, scmax=0, nacol="grey80") {
    if (length(color)==1) color <- c(color, "white", color)
    if (length(color)==2) color <- c(color[1], "white", color[2])
    if (scmax==0) scmax <- max(abs(x), na.rm=T)
    pos <- which(abs(x) > scmax)
    if (length(pos)>0) x[pos] <- scmax*sign(x[pos])
    x <- abs(x/scmax)^gama*sign(x)
    color <- t(col2rgb(color))
    col <- sapply(x, function(x, color) {
        colSums(color*c(abs(x)*(x<0), 1-abs(x), x*(x>0)))
    }, color=color/255)
    pos <- which(colSums(is.na(col))>0)
    col[is.na(col)] <- 0
    col <- apply(col, 2, function(x, alpha) rgb(x[1], x[2], x[3], alpha=alpha), alpha=alpha)
    col[pos] <- nacol
    return(col)
}

#' Plot heatmap
#'
#' This function produce a heatmap plot from a numerical matrix
#'
#' @param x Numerical matrix
#' @param color Two character strings vector describing the colors for the heatmap
#' @param gama Number, indicating the exponential transformation for the color scale
#' @param cex Number indicating the magnification factor for the labels
#' @param grid Logical, whether a grid should be ploted
#' @param scale Number between 0 and .9 indicating the proportion of vertical space used to draw the color scale
#' @param scmax Optional number indicating the maximum value to be allowed for the heatmap
#' @param box Logical, whether to draw a box around the plot
#' @param ... Additional parameters to pass to the plot function
#' @return Nothing, a heatmap is produced in the default output device
#' @export

plothm <- function(x, color=c("royalblue","firebrick2"), gama=1, cex=1, grid=T, scale=F, scmax=0, box=TRUE, ...) {
    if (scale>0) {
        if (scale==1) ff <- 6/(nrow(x)+5)
        else ff <- scale
        pari <- par("mai")
        layout(matrix(1:2, 2, 1), h=c(1-ff, ff))
        if (round(sum(pari-c(1.02, .82, .82, .42)), 2)==0) pari <- c(.2, .2, 1.2, 1.2)
        par(mai=pari)
        plothm.(x, color=color, gama=gama, scmax=scmax, box=box, ...)
        axis(4, nrow(x):1, rownames(x), tick=F, line=0, las=2, adj=0, cex.axis=cex)
        axis(3, 1:ncol(x), colnames(x), tick=F, line=0, las=2, adj=0, cex.axis=cex)
        ra <- seq(-1, 1, length=100)
        coli <- colorScale(x=ra, color=color, gama=gama, scmax=scmax)
        par(mai=pari*c(0, 1, 0, 3)+c(.5, 0, .1, 0))
        image(1:length(ra), 1, matrix(1:length(ra), length(ra), 1), col = coli, ylab = "", xlab = "", axes = F)
        if (scmax==0) scmax <- max(abs(x), na.rm=T)
        axis(1, seq(1, length(ra), length=5), round(seq(-scmax, scmax, length=5), 1), cex.axis=cex)
    }
    else plothm.(x=x, color=color, gama=gama, grid=grid, scmax=scmax, box=box, ...)
}