<p>I am trying to plot a heatmap of read density around a feature of interest (TSSs) very common in genomics papers. something like this (B):</p>

<p><img src='http://ars.els-cdn.com/content/image/1-s2.0-S0092867410003181-gr3.jpg' alt='enter image description here' /></p>

<p>However, I am struggling a bit in getting to look "right". A bit of  background:
I have mapped ChIP-seq reads for pol2 and calculate the coverage, per nucleotide, using <a href='https://github.com/arq5x/bedtools2'>bedtools</a>. </p>

<pre><code>coverageBed -d -abam $bamFile -b $TSSs &gt; $coverage.bed
# output:
chr1    67108226    67110226    uc001dct.3    16    +    1    10
chr1    67108226    67110226    uc001dct.3    16    +    2    10
chr1    67108226    67110226    uc001dct.3    16    +    3    10
chr1    67108226    67110226    uc001dct.3    16    +    4    10
chr1    67108226    67110226    uc001dct.3    16    +    5    8
chr1    67108226    67110226    uc001dct.3    16    +    6    8
chr1    67108226    67110226    uc001dct.3    16    +    7    8
chr1    67108226    67110226    uc001dct.3    16    +    8    8
chr1    67108226    67110226    uc001dct.3    16    +    9    8
chr1    67108226    67110226    uc001dct.3    16    +    10    8
</code></pre>

<p>Then in R, the genomic position, in column 7, is converted to relative position to the TSS and read counts normalized to the library size. This is converted to a numeric matrix with each row being a TSS and each column the relative nucleotide position. For the plotting the matrix is ordered number of reads per TSS, and the values logged. This is the outcome:</p>

<pre><code>heatmap(cov.mlog, Rowv=NA, Colv=NA, scale="column", labCol = FALSE, labRow = FALSE, col=brewer.pal(9, "Greens"), margins = c(5, 5))
</code></pre>

<p><img src='https://dl.dropboxusercontent.com/u/9519660/testheatmap1.jpg' alt='enter image description here' /></p>

<p>Although the average coverage plot looks as one would expect for pol2, the heatmap is, well, not great. My question is what am I doing wrong and how to improve it? In <a href='http://www.sciencedirect.com/science/article/pii/S0092867410003181'>this</a> and <a href='http://genomebiology.com/2012/13/10/r85#sec4'>this</a> paper, the coverage is calculate for a bin of nucleotides and not per nucleotide as I am doing. Would that improve visualization? How to do it in <a href='https://github.com/arq5x/bedtools2'>bedtools</a>? Should the sorting of the matrix be done differently?</p>

<p>I am aware that could use ngsplot for this, but I am trying to avoid it because my implementation would fit better with my other analysis.</p>

<p>Thank you!</p>

<p>Hi,</p>

<p>you may find <a href='https://github.com/fidelram/deepTools'>deeptools</a> useful. </p>

<p>Example plots can be found here: <a href='http://f1000.com/posters/browse/summary/1094053'>http://f1000.com/posters/browse/summary/1094053</a></p>

<p>Furthermore, deeptools can create bigwig files that are normalized or ratios. As a bonus a Galaxy wrapper is available in the <a href='http://testtoolshed.g2.bx.psu.edu/view/bgruening/deeptools'>ToolShed</a>. </p>
