Time and memory efficient processing of many intervals for a bigWig file

Question

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.

Answer 1

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