title: GenomicRanges add coverage
body: <p>I'm working on RNA seq data and trying to plot average coverage profiles by genotype, similarly to what is done here</p>
  <p><strong>RNA seq Coverage per genotype (Source: pickrell et al, Nature, 2010)</strong></p>
  <p><a href="https://i.stack.imgur.com/KHDif.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/KHDif.png" alt="enter image description here" /></a></p>
  <p>To do this plot, I have bigwig files from 100 individuals, that contain coverage information from RNA-seq data (in a specific region) and that I read in R, as GenomicRanges objects.</p>
  <p>this gives me GRanges objects such as those obtained in the following toy example :</p>
  <blockquote>
  <p>gr1=GRanges(seqname=1,range=IRanges(start=c(1,5,10,15,30,55), end=c(4,9,14,29,39,60)))</p>
  <p>gr1$cov=c(3,1,8,6,2,10)</p>
  <p>gr2=GRanges(seqname=1,range=IRanges(start=c(3,20,24), end=c(7,23,26)))</p>
  <p>gr2$cov=c(3,5,3)</p>
  <p>start=unique(sort(c(ranges(gr1)@start,ranges(gr2)@start)))</p>
  <p>gr1</p>
  </blockquote>
  <pre><code>GRanges object with 6 ranges and 1 metadata column:
  seqnames    ranges strand |       cov
     &lt;Rle&gt; &lt;IRanges&gt;  &lt;Rle&gt; | &lt;numeric&gt;
         1  [ 1,  4]      * |         3
         1  [ 5,  9]      * |         1
         1  [10, 14]      * |         8
         1  [15, 29]      * |         6
         1  [30, 39]      * |         2
         1  [55, 60]      * |        10 
          -------
   seqinfo: 1 sequence from an unspecified genome; no seqlengths
  </code></pre>
  <blockquote>
  <p>gr2</p>
  </blockquote>
  <pre><code>GRanges object with 3 ranges and 1 metadata column:
  seqnames    ranges strand |       cov
     &lt;Rle&gt; &lt;IRanges&gt;  &lt;Rle&gt; | &lt;numeric&gt;
         1  [ 3,  7]      * |         3
         1  [20, 23]      * |         5
         1  [24, 26]      * |         3
         -------
   seqinfo: 1 sequence from an unspecified genome; no seqlengths
  </code></pre>
  <p>The problem is that I have these per individual (gr1 and gr2 would be 2 different individuals), and I'd like to combine them to create a genomic ranges object that gives me the total coverage at each position across individuals, 1 and 2
  that would look as follows :</p>
  <blockquote>
  <p>gr3</p>
  </blockquote>
  <pre><code>GRanges object with 6 ranges and 1 metadata column:
  seqnames    ranges strand |       cov
     &lt;Rle&gt; &lt;IRanges&gt;  &lt;Rle&gt; | &lt;numeric&gt;
         1  [ 1,  2]      * |         3
         1  [ 3,  4]      * |         6 (=3+3)
         1  [ 5,  7]      * |         4 (=1+3)
         1  [ 8,  9]      * |         1
         1  [10, 14]      * |         8
         1  [15, 19]      * |         6
         1  [20, 23]      * |         11 (=6+5)
         1  [24, 26]      * |         9 (=6+3)
         1  [27, 29]      * |         6
         1  [30, 39]      * |         2
         1  [55, 60]      * |        10 
  </code></pre>
  <p>Does anyone know a simple way to do this ? or am I doomed ?</p>
  <p>Thanks for your answers.</p>
  <p>PS:
  my data isn't stranded, but if you have it for stranded data it's even better.</p>
  <p>PPS: Ideally, I would also like to be able to to a multiplication, or apply any function with two arguments x and y, instead of simply adding the coverage.</p>