Skip to content

Instantly share code, notes, and snippets.

@jokergoo
Last active December 19, 2015 06:09
Show Gist options
  • Save jokergoo/5909943 to your computer and use it in GitHub Desktop.
Save jokergoo/5909943 to your computer and use it in GitHub Desktop.
use heatmap to represent density distribution
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)
@jokergoo
Copy link
Author

jokergoo commented Jul 2, 2013

density

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment