Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T22:21:28.360Z
Answers.Text: Well, this looks very promising! I've only come across bx-python before for accessing BigWig files before.
Answers.UserDisplayName: Ian Sudbery
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-14T06:44:59.860Z
Answers.Text: It should also be a bit faster than bx-python (I wrote it largely to speed up deepTools).
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-14T13:06:26.420Z
Answers.Text: 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.
Answers.UserDisplayName: Ian Sudbery
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-14T13:08:59.647Z
Answers.Text: @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.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-15T09:56:29.397Z
Answers.Text: 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.
Answers.UserDisplayName: Ian Sudbery
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-15T11:12:31.473Z
Answers.Text: The format is optimized for using zoom levels, which you're not using in this case.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T20:43:29.640Z
Answers.LastActivityDate: 2017-06-13T20:52:36.063Z
Answers.LastEditDate: 2017-06-13T20:52:36.063Z
Answers.UserDisplayName: Devon Ryan
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-14T10:17:00.347Z
Comments.Text: 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.
Comments.UserDisplayName: Konrad Rudolph
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-14T12:05:30.163Z
Comments.Text: 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.
Comments.UserDisplayName: Ian Sudbery
Tags: <bigwig>
Title: Time and memory efficient processing of many intervals for a bigWig file
ThreadText: # 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](http://genesdev.cshlp.org/content/29/1/63.abstract).
  
  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](http://genome.ucsc.edu/cgi-bin/hgFileUi?db=hg19&g=wgEncodeMapability). 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:
  
  ```r
  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**: @username_0: 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.  
  
  ---
  
  