plot-gene-histogram.R

## Plot heatmap in the style of Fig 1B-2:
## dx.doi.org/10.1371/journal.pgen.1006224.g001
library(ggplot2)

## Mock data of 50 rows, 10 columns and 2 strands.
value <- function(mean_ = 0, length_ = 10) {
    dnorm(seq(-1, 3, length = length_), mean = mean_, sd = 0.2)
}
noise <- function(length_ = 10) {
    rnorm(length_, mean = 0, sd = 0.1)
}
normalize <- function(x) {
    (x - min(x)) / (max(x) - min(x))
}
genes <- data.frame(
    enum = rep(1:50, each = 10),
    pos = rep.int(seq(-1000, 3000, length = 10), 50),
    plus = c(replicate(50, value() + noise())),
    minus = c(sapply(seq(0, 3, length = 50), value) + replicate(50, noise()))
)
## ^ Normalize counts to [0, 1] to honor alpha channel limits.

## Per http://stackoverflow.com/a/14496538
ggplot(genes, aes(pos, enum)) + theme_bw() +
    coord_cartesian(expand = FALSE) +
    scale_y_reverse() +
    geom_tile(alpha = normalize(genes$plus), fill = "red") +
    geom_tile(alpha = normalize(genes$minus), fill = "blue")