<p>Assuming you can make a BED file of introns, you can then use the <code>pyBigWig</code> module in python:</p>

<pre><code>import pyBigWig
import numpy

bw = pyBigWig.open("some file.bw")
bed = open("introns.bed")
for line in bed:
    cols = line.strip().split("\t")
    vals = bw.values(cols[0], int(cols[1]), int(cols[2]), numpy=True)
    effLen = cols[2] - cols[1] - (vals &lt; 1.0).sum()
    # Do something with this.
</code></pre>

<p>One can make that parallel (with the deeptools API) should it prove too slow. This will only use an appreciable amount of memory for very large introns, since there you need to store the value at each base. One could make this much more memory efficient by using <code>bw.intervals()</code> instead, but that'd require quite a bit more code (you'd need to determine the number of bases overlap per interval and then assume that either each base you want to query has an overlapping interval or, if not, keep track of that fact).</p>

CC BY-SA 3.0
2017-06-13T22:21:28.360Z
Well, this looks very promising! I've only come across bx-python before for accessing BigWig files before.
Ian Sudbery
CC BY-SA 3.0
2017-06-14T06:44:59.860Z
It should also be a bit faster than bx-python (I wrote it largely to speed up deepTools).
Devon Ryan
CC BY-SA 3.0
2017-06-14T13:06:26.420Z
This solution takes about 2 hours to process the file on our system in case you are interested. Although I wasn't able to use the numpy version. Uses very little memory as far as I can tell.
Ian Sudbery
CC BY-SA 3.0
2017-06-14T13:08:59.647Z
@IanSudbery: Thanks for the update. I don't expect the numpy version to be appreciably faster. If you need this faster then I can make a multithreaded version using the deepTools API.
Devon Ryan
CC BY-SA 3.0
2017-06-15T09:56:29.397Z
Out of interest, why is BigWig access so slow. If this were a BAM file or a tabix indexed bedgraph, the process would have taken seconds with psyam or Rsamtools or similar.
Ian Sudbery
CC BY-SA 3.0
2017-06-15T11:12:31.473Z
The format is optimized for using zoom levels, which you're not using in this case.
Devon Ryan
CC BY-SA 3.0
2017-06-13T20:43:29.640Z
2017-06-13T20:52:36.063Z
2017-06-13T20:52:36.063Z
Devon Ryan
<p>I am trying to calculate the mappability adjusted length of introns as described by <a href="http://genesdev.cshlp.org/content/29/1/63.abstract" rel="nofollow noreferrer">Boutz et al</a>.</p>

<p>Briefly, for each intron I wish to calculate the length minus the number of bases that are non-uniquely mappable. Mappability tracks can be downloaded as bigWigs from <a href="http://genome.ucsc.edu/cgi-bin/hgFileUi?db=hg19&amp;g=wgEncodeMapability" rel="nofollow noreferrer">here</a>. The the score at each base is 1/number of mapping positions, so 1 indicates a read can be uniquely mapped at this location. </p>

<p>I am struggling to do this computation in either a reasonable amount of time, or a reasonable amount of memory.</p>

<p>First I tried importing the whole bigwig:</p>

<pre><code>library(rtracklayer)
mappability &lt;- import(BigWigFile("wgEncodeCrgMapabilityAlign50mer.bigWig"),
                      as="NumericList",
                      selection=BigWigSelection(intron_ranges))
</code></pre>

<p>where <code>intron_ranges</code> is a <code>GRanges</code> object with the introns (about 800,000 of them).</p>

<p>This used <strong>way</strong> too much memory, and soon caused by machine to fall over.</p>

<p>Second I tried processing one intron at a time:</p>

<pre class="lang-r prettyprint-override"><code>mappability_file = BigWigFile("wgEncodeCrgMapabilityAlign50mer.bigWig")
effective_length &lt;- function (gr) {
  intron_selection &lt;- BigWigSelection(gr)
  scores &lt;- import(mappability_file, as="NumericList", selection=intron_selection)
  non_unique &lt;- sum(scores[[1]] &lt; 1.0)
  eff_len = width(gr)[1] - non_unique
  return(eff_len)
}    

mappability &lt;- sapply(intron_ranges, effective_length)
</code></pre>

<p>This has been running for hours and shows no sign of finishing. </p>

<p>Is there a way to do this that uses less than, say 4GB of RAM, but finishes in say less than 30 minutes? It feels like there should be. </p>

<p>I'm not tied to R, just what I've tried so far. Happy with answers using binary packages, common shell tools, python or R. </p>

CC BY-SA 3.0
2017-06-14T10:17:00.347Z
4 GiB is *not* a lot of memory, especially when working with genomic data. In fact, 800.000 `GRanges` isn’t all that much (I’m working with repeat elements, I have *a lot* more). … What I’m trying to say: time to upgrade your hardware.
Konrad Rudolph
CC BY-SA 3.0
2017-06-14T12:05:30.163Z
I have 32GB on my machine. Its not the Granges that takes the space, its the passing them through import.bw. I siad 4GB as I don't want to have to close everything else down to do this sort of thing (processing BigWigs must be a pretty common task). I have access to a machine with 256GB on the local cluster so it would be possible to use that, but it seems unnecessary. Genome aligners that store indexes in memory use less that than. Even STAR which store all sorts fits in 24GB.
Ian Sudbery
<bigwig>
Time and memory efficient processing of many intervals for a bigWig file
# Time and memory efficient processing of many intervals for a bigWig file

(Asked by: username_0 on 2017-06-13)

I am trying to calculate the mappability adjusted length of introns as described by Boutz et al.
Briefly, for each intron I wish to calculate the length minus the number of bases that are non-uniquely mappable. Mappability tracks can be downloaded as bigWigs from here. The the score at each base is 1/number of mapping positions, so 1 indicates a read can be uniquely mapped at this location. 
I am struggling to do this computation in either a reasonable amount of time, or a reasonable amount of memory.
First I tried importing the whole bigwig:
library(rtracklayer)
mappability <- import(BigWigFile("wgEncodeCrgMapabilityAlign50mer.bigWig"),
                      as="NumericList",
                      selection=BigWigSelection(intron_ranges))

where intron_ranges is a GRanges object with the introns (about 800,000 of them).
This used way too much memory, and soon caused by machine to fall over.
Second I tried processing one intron at a time:
mappability_file = BigWigFile("wgEncodeCrgMapabilityAlign50mer.bigWig")
effective_length <- function (gr) {
  intron_selection <- BigWigSelection(gr)
  scores <- import(mappability_file, as="NumericList", selection=intron_selection)
  non_unique <- sum(scores[[1]] < 1.0)
  eff_len = width(gr)[1] - non_unique
  return(eff_len)
}    

mappability <- sapply(intron_ranges, effective_length)

This has been running for hours and shows no sign of finishing. 
Is there a way to do this that uses less than, say 4GB of RAM, but finishes in say less than 30 minutes? It feels like there should be. 
I'm not tied to R, just what I've tried so far. Happy with answers using binary packages, common shell tools, python or R.

**Tags:** bigwig
> **Comments:**
> **username_1**: 4 GiB is *not* a lot of memory, especially when working with genomic data. In fact, 800.000 `GRanges` isn’t all that much (I’m working with repeat elements, I have *a lot* more). … What I’m trying to say: time to upgrade your hardware.  
> **username_0**: I have 32GB on my machine. Its not the Granges that takes the space, its the passing them through import.bw. I siad 4GB as I don't want to have to close everything else down to do this sort of thing (processing BigWigs must be a pretty common task). I have access to a machine with 256GB on the local cluster so it would be possible to use that, but it seems unnecessary. Genome aligners that store indexes in memory use less that than. Even STAR which store all sorts fits in 24GB.  

---

## 1 answers

### Answer by username_2 on 2017-06-13. Score: 2 | Accepted answer

Assuming you can make a BED file of introns, you can then use the pyBigWig module in python:
import pyBigWig
import numpy

bw = pyBigWig.open("some file.bw")
bed = open("introns.bed")
for line in bed:
    cols = line.strip().split("\t")
    vals = bw.values(cols[0], int(cols[1]), int(cols[2]), numpy=True)
    effLen = cols[2] - cols[1] - (vals < 1.0).sum()
    # Do something with this.

One can make that parallel (with the deeptools API) should it prove too slow. This will only use an appreciable amount of memory for very large introns, since there you need to store the value at each base. One could make this much more memory efficient by using bw.intervals() instead, but that'd require quite a bit more code (you'd need to determine the number of bases overlap per interval and then assume that either each base you want to query has an overlapping interval or, if not, keep track of that fact).

> **Comments:**
> **username_0**: Well, this looks very promising! I've only come across bx-python before for accessing BigWig files before.  
> **username_2**: It should also be a bit faster than bx-python (I wrote it largely to speed up deepTools).  
> **username_0**: This solution takes about 2 hours to process the file on our system in case you are interested. Although I wasn't able to use the numpy version. Uses very little memory as far as I can tell.  
> **username_2**: @IanSudbery: Thanks for the update. I don't expect the numpy version to be appreciably faster. If you need this faster then I can make a multithreaded version using the deepTools API.  
> **username_0**: Out of interest, why is BigWig access so slow. If this were a BAM file or a tabix indexed bedgraph, the process would have taken seconds with psyam or Rsamtools or similar.  
> **username_2**: The format is optimized for using zoom levels, which you're not using in this case.  

---

