How to extract bigWig signal for a given bed file?

Question

[Bioradical]

I'm looking for a program to extract signal values from a bigWig file for a given set of coordinates.

Python, R or command line is fine (my programming is rudimentary at best). Thank you!

Answer 1

[Devon Ryan]

In python:

import pyBigWig
bw = pyBigWig.open("some_file.bw")
for line in open("foo.bed"):
    cols = line.strip().split()
    vals = bw.values(cols[0], int(cols[1]), int(cols[2]))
    # Do something with the values...
bw.close()

[kashiff007]

cool library... does this also takes care about negative and positive strands from bed file?

[Devon Ryan]

There's nothing special about the minus strand, it has the same values in a bigWig file.

[kashiff007]

I think it will be helpful if one is plotting the nucleosome signal around TSS or TTS. As some genes has reverse orientation (negative strands in gff file), there data needs to be reversed accordingly before plotting. Correct me if I am wrong.

[Devon Ryan]

Tools like computeMatrix in deepTools handle such things automatically (deepTools uses the pyBigWig library under the hood).

[kashiff007]

Ohh... thanks I have used it, but completely forgotten.

[John]

Smooth as hell :) Great work Devon!

[GenoMax]

Are you done with the thesis?

[Devon Ryan]

I was going to yell that just now when I saw him walk into a talk :)

[GenoMax]

Looks like @John wants to stay a student .. forever.

Not smooth, @John :)

[John]

Hehehe, hasn't there been enough bioinformatician-abuse for one day? :P

For the record, my answer would have been "uhhh..."

[Devon Ryan]

You're a grad student, there's never enough grad student abuse :)

[mccormack]

Hi Devon, Thank you very much for this script. I added three lines where you have the comment, 'Do something with the values':

ATG =  cols[3]           # get the identifier from the bed file 
span = sum(vals)         # sum all the values in vals
print(ATG,"\t",span,"\n")

This prints out a list of ATG numbers and a number. For example,

AT5G27350 6375.189992427826

However, some ATG numbers have 'nan' instead of a number, such as :

AT5G64667 nan

I was wondering if the reason why 'nan' is printed is because that in the list there may be 1 or 2 nan's, so the sum for the whole list is 'nan' ?
When I put 'nan' in a list I created and try to sum all value in this list with sum(mylist), I get an error, so I don't know how 'nan' is working in the bw.values() function.
If there is just a single 'nan' in the vals, then does that mean a sum of the vals would be 'nan' ? If this is true, is there a way to replace the 'nan' with a 0 ?

[Devon Ryan]

You'll want np.nansum().

[mccormack]

Great ! Thank you very much. Works very well. Numbers for all the ATGs (genes).

[hani.kim127]

Hi, I have been struggling to find a way to extract signal from a bw file. And this code is perfect for what I want to do. However, I do not know python. I would like to do this on Linux terminal.

Could you share the code that replicates the above in linux?

Thanks

[Devon Ryan]

Type that into a file, pip install --user pyBigWig and then execute the file you saved.

Answer 2

[dariober]

Altrenatively, use bigWigToBedGraph from UCSC genome browser at http://hgdownload.soe.ucsc.edu/admin/exe/, just pick the compiled version that suites you. Usage:

bigWigToBedGraph - Convert from bigWig to bedGraph format.
usage:
   bigWigToBedGraph in.bigWig out.bedGraph
options:
   -chrom=chr1 - if set restrict output to given chromosome
   -start=N - if set, restrict output to only that over start
   -end=N - if set, restict output to only that under end
   -udcDir=/dir/to/cache - place to put cache for remote bigBed/bigWigs

[Bioradical]

In bedgraph format, the fourth column is == to the signal?

[Devon Ryan]

Yes, it's whatever signal is stored in the bigWig file.

Answer 3

[chechu086]

You also can use the BigWigAverageOverBed program from UCSC. The usage: bigWigAverageOverBed in.bw in.bed out.tab