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:

mappability <- import(BigWigFile("wgEncodeCrgMapabilityAlign50mer.bigWig"),

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

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.

  • 1
    $\begingroup$ 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. $\endgroup$ – Konrad Rudolph Jun 14 '17 at 10:17
  • $\begingroup$ 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. $\endgroup$ – Ian Sudbery Jun 14 '17 at 12:05

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).

  • $\begingroup$ Well, this looks very promising! I've only come across bx-python before for accessing BigWig files before. $\endgroup$ – Ian Sudbery Jun 13 '17 at 22:21
  • $\begingroup$ It should also be a bit faster than bx-python (I wrote it largely to speed up deepTools). $\endgroup$ – Devon Ryan Jun 14 '17 at 6:44
  • 1
    $\begingroup$ 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. $\endgroup$ – Ian Sudbery Jun 14 '17 at 13:06
  • $\begingroup$ @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. $\endgroup$ – Devon Ryan Jun 14 '17 at 13:08
  • 1
    $\begingroup$ The format is optimized for using zoom levels, which you're not using in this case. $\endgroup$ – Devon Ryan Jun 15 '17 at 11:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.