use heatmap to represent density distribution

density.heatplot.R

require(RColorBrewer)
 
 
# since gene expression data always contains a lot of samples.
# If one want to see distributions of samples, using boxplot would be not so clear.
# The function use continuous colors to represent density distributions of expression
# values in samples and can make a better illustration of the data.
heatplot = function(x, col=rev(brewer.pal(10, "Spectral")), draw.quantiles = TRUE, align = TRUE, each = FALSE, ...) {
    if(is.vector(x) && class(x) != "list") {
        x = as.matrix(x)
    }
    if(is.matrix(x)) {
        n = dim(x)[2]
        # if different styles of colors are used, it should be formatted as a list
        # because the number of sections of colors may be different
 
        dx = apply(x, 2, function(x) density(x)$x)  # data value
        dy = apply(x, 2, function(x) density(x)$y)  # density value
        quantile.values = apply(x, 2, quantile)
        mean.values = apply(x, 2, mean)
    }
    if(is.list(x)) {
        n = length(x)
 
        dx = sapply(x, function(x) density(x)$x)  # data value
        dy = sapply(x, function(x) density(x)$y)  # density value
        quantile.values = sapply(x, quantile)
        mean.values = sapply(x, mean)
    }
 
    if(!is.list(col)) {
        col = rep(list(col), n)
    }
    if(is.list(col) && length(col) != n) {
        stop("Since 'col' is specified as a list, it should has the same length as numbers of columns in 'x'.")
    }
    if(!all(sapply(col, length) > 1)) {
        stop("Length of any color vector should contain at least two colors.")
    }
    
    if(! each) {
        min.density = min(as.vector(dy))
        max.density = max(as.vector(dy))
        range.density = max.density - min.density
        dy = (dy - min.density) / range.density
    }

    min.value = min(as.vector(dx))
    max.value = max(as.vector(dx))
    range.value = max.value - min.value
 
    plot(c(0, n+1), c(min.value, max.value), type = "n", axes=FALSE, ann=FALSE, ...)
    for(j in 1:n) {

        if(each) {
            min.density = min(dy[, j])
            max.density = max(dy[, j])
            range.density = max.density - min.density
            dy[, j] = (dy[, j] - min.density) / range.density
        }

        for(i in 2:length(dy[, j])) {
            color = color.pal(dy[i, j], col=col[[j]], breaks=seq(0, 1, length.out=length(col[[j]])))
            rect(j-0.5, dx[i-1, j], j+0.5, dx[i, j], col=color, border=color)
        }
 
        if(align) {
            color = color.pal(min.density, col=col[[j]], breaks=seq(0, 1, length.out=length(col[[j]])))
            rect(j-0.5, min(dx[, j]), j+0.5, min.value, col=color, border=color)
            rect(j-0.5, max(dx[, j]), j+0.5, max.value, col=color, border=color)
        }
    }
    #axis(side = 2)
    if(draw.quantiles) {
        for(i in 1:dim(quantile.values)[1]) {
            lines(1:n, quantile.values[i, ], col="black", lwd=1)
        }
        lines(1:n, mean.values, col = "black", lwd = 1)
        text(rep(n+0.6, dim(quantile.values)[1]), quantile.values[, n], rownames(quantile.values), cex=0.8, adj=c(0, 0.5))
        text(n+0.6, mean.values[n], "mean", cex=0.8, adj=c(0, 0.5))
    }
}

jitplot = function(x, alpha = 0.05) {

    if(is.matrix(x)) {
        n = dim(x)[2]
    }

    min.value = min(as.vector(x))
    max.value = max(as.vector(x))
    range.value = max.value - min.value

    plot(c(0, n+1), c(min.value, max.value), type = "n", axes=FALSE, ann=FALSE)
    for(j in 1:n) {
        k = length(x[, j])
        points((runif(k)-0.5)*0.8+j, x[, j], col = rgb(0, 0, 0, alpha), pch=16)
    }

}

x = NULL

for(i in 1:30) {
    mean = runif(1)*5
    sd = runif(1)+5
    x = cbind(x, rnorm(1000, mean, sd))
}

par(mfrow = c(2, 2), mar = c(1, 1, 1, 1))
heatplot(x, each = TRUE, draw.quantiles = FALSE)
jitplot(x)
vioplot(x[,1],x[,2],x[,3],x[,4],x[,5],x[,6],x[,7],x[,8],x[,9],x[,10],
        x[,11],x[,12],x[,13],x[,14],x[,15],x[,16],x[,17],x[,18],x[,19],x[,20],
        x[,21],x[,22],x[,23],x[,24],x[,25],x[,26],x[,27],x[,28],x[,29],x[,30])
boxplot(x)