07-plot_heatmap.Rmd

# Plot heatmap

```{r Plot heatmap of plus strand cropped to median gene size, cache = TRUE}
load("reads_h3u.RData")
load("counts.RData")

genes_plot <- granges(counts$full)
## Subset to positive strand.
genes_plot <- genes_plot[strand(genes_plot) == "+"]
## Crop to median size to reduce plot width.
summary(width(genes_plot))
genes_plot_long <- genes_plot
width_median <- median(width(genes_plot_long))
width(genes_plot) <- pmin(width(genes_plot_long), as.integer(width_median))
stopifnot(all(start(genes_plot) == start(genes_plot_long)))
save(genes_plot, file = "genes_plot.RData")
summary(width(genes_plot))

## Create tiles to generate count matrix.
window_size <- 500
windows <- tile(genes_plot, width = window_size)
## Create count matrix.
scoreOverlaps <- function(gr, reads) {
    ol <- findOverlaps(gr, reads)
    mcols(ol)$score <- score(reads[to(ol)])
    df <- as(ol, "DataFrame")
    df <- aggregate(score ~ queryHits, df, sum)
    score <- vector("integer", length(gr))
    score[df$queryHits] <- df$score
    score
}
counts_list <- relist(scoreOverlaps(unlist(windows), reads_h3u), windows)
mat <- as(counts_list, "matrix")
## Sort by window counts.
mat[is.na(mat)] <- 0
mat <- mat[order(rowSums(mat)), ]
tail(mat)[, 1:6]
## Plot.
image(x = 1:ncol(mat) * window_size,
      y = 1:nrow(mat),
      z = t(log2(mat + 1)),
      main = "Human 3 untreated (+ strand)",
      xlab = "Distance from TSS (bp)",
      ylab = "Gene")
```