Question: How to extract bigWig signal for a given bed file?
1
gravatar for Bioradical
2.2 years ago by
Bioradical50
United States
Bioradical50 wrote:

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!

bigwig • 2.2k views
ADD COMMENTlink modified 2.2 years ago by dariober9.7k • written 2.2 years ago by Bioradical50
5
gravatar for Devon Ryan
2.2 years ago by
Devon Ryan86k
Freiburg, Germany
Devon Ryan86k wrote:

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()
ADD COMMENTlink written 2.2 years ago by Devon Ryan86k

Smooth as hell :) Great work Devon!

ADD REPLYlink written 2.2 years ago by John12k
1

Are you done with the thesis?

ADD REPLYlink modified 2.2 years ago • written 2.2 years ago by genomax59k
1

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

ADD REPLYlink written 2.2 years ago by Devon Ryan86k
1

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

Not smooth, @John :)

ADD REPLYlink modified 2.2 years ago • written 2.2 years ago by genomax59k

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

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

ADD REPLYlink modified 2.2 years ago • written 2.2 years ago by John12k
1

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

ADD REPLYlink written 2.2 years ago by Devon Ryan86k

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

ADD REPLYlink written 9 months ago by kashifalikhan00740

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

ADD REPLYlink written 9 months ago by Devon Ryan86k

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.

ADD REPLYlink written 9 months ago by kashifalikhan00740

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

ADD REPLYlink written 9 months ago by Devon Ryan86k

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

ADD REPLYlink written 9 months ago by kashifalikhan00740

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 ?

ADD REPLYlink modified 12 days ago • written 12 days ago by mccormack20

You'll want np.nansum().

ADD REPLYlink written 12 days ago by Devon Ryan86k

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

ADD REPLYlink written 12 days ago by mccormack20
2
gravatar for dariober
2.2 years ago by
dariober9.7k
Glasgow - UK
dariober9.7k wrote:

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
ADD COMMENTlink written 2.2 years ago by dariober9.7k

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

ADD REPLYlink written 2.2 years ago by Bioradical50

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

ADD REPLYlink written 2.2 years ago by Devon Ryan86k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1157 users visited in the last hour