This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Identifying The 'Summit' Coordinate Within Many Coordinates In A Wig File?

I have a WIG file containing read coverage data from ChIP-seq peak calling analysis.

For a given set of coordinates is there a simple fast way of determining the 'summit' coordinate, i.e the coordinate with the highest read pileup.

My own method would be slow as i would run through the (big) WIG file, extract the region of interest and find the coordinate with the highest score. But surely there must be a quicker way of doing this? If not fair enough.

The peak caller does report the summit from its own data, but i want to check out different regions.

Thanks.

wiggle
• 6,756 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by Ian Ian

2 answers

Convert your Wig file to BigWig with the UCSC executable wigToBigWig. This provides an indexed Wig file you can query to return regions of interest, without having to loop over the entire file each time.

bx-python provides a python BigWig interface that will help with finding your summit coordinates:

from bx.bbi.bigwig_file import BigWigFile

f = open("your_file.bw")
bw = BigWigFile(file=f)
vals = bw.get("chr3", 1, 1000)
# vals is a list of start, end, value
max_start, max_end, max_val = max(vals, key=lambda x: x[-1])

You'll want your logic to be smarter to handle multiple summits of the same size but hopefully that helps to get started.

• 0 views • link
written by Brad Chapman Brad Chapman

Looks interesting! Shame its in python though as i am a Perl'ist :)

• 0 views • link
written by Ian Ian

convert your wig file to bigwig using wig2bigwig and use bigWigSummary to extract the information in a given region.

   bigWigSummary file.bigWig chrom start end dataPoints

Get summary data from bigWig for indicated region, broken into
dataPoints equal parts.  (Use dataPoints=1 for simple summary.)

options:
   -type=X where X is one of:
         mean - average value in region (default)
         min - minimum value in region
         max - maximum value in region
         std - standard deviation in region
         coverage - % of region that is covered
• 270 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thanks. I was aware of this, but it only returns the value, not the coordinate.

• 0 views • link
written by Ian Ian

Log in to answer this question.

More posts like this