bigWig to bed for regions above/below threshold

Question

[skanterakis]

I have a bigWig file with the following profile:

bigWigInfo uniqueness.bw
version: 4
isCompressed: yes
isSwapped: 0
primaryDataSize: 919,783,047
primaryIndexSize: 97,088,256
zoomLevels: 10
chromCount: 25
basesCovered: 3,095,693,183
mean: 0.828284
min: 0.000000
max: 1.000000
std: 0.361063

I would like to create a bed for all regions with value < 1.

More generally, I'd like to be able to threshold a bigWig on a value and get all regions above or below that value in bed format.

Thanks a lot!

[Pierre Lindenbaum]

Hello skanterakis!

It appears that your post has been cross-posted to another site: https://bioinformatics.stackexchange.com/questions/2800/

This is typically not recommended as it runs the risk of annoying people in both communities.

[oriolebaltimore]

Dear Alex and Sean, could you please help understand what is the threshold value mean in either ATAC-Seq or ChIP-Seq in wig file obtained after bigWigToWig conversion. Particularly, in bedgraph format what is the 4th column in dataValue and what are the norms for selecting the threshold. Thanks Adrian

[Devon Ryan]

There is no fixed meaning to that value.

Answer 1

[Devon Ryan]

Copy-pasting from where you cross-posted:

You can do that with a bit of python:

#!/usr/bin/env python
import pyBigWig
threshold = 10  # Change me
bw = pyBigWig.open("file.bw")  # Change me
of = open("regions.bed", "w")  # Change me

for chrom, len in bw.chroms().items():
    intervals = bw.intervals(chrom)
    for interval in interals:
        if abs(interval[2]) > threshold:
            of.write("{}\t{}\t{}\n".format(chrom, interval[0], interval[1]))
bw.close()
of.close()

You'll need to install pyBigWig (it's available via pip or conda install) and change the lines with Change me in them.

Answer 2

[Pierre Lindenbaum]

I would like to create a bed for all regions with value < 1.

convert the wig back to bed: Any Method Of Converting Bigwig File Format Into Bed Format? , filter with awk ( awk '($4< 1.0)' ) , convert back to bigwig if needed with bedgraphtobigwig

Answer 3

[Sean Davis]

1. Use import() from the rtracklayer bioconductor package to read the bigwig file.
2. Use slice() from the GenomicRanges bioconductor package on the results from #1.
3. Use export() from the rtracklayer bioconductor package on the results from #2 to write the BED file.

Answer 4

[Alex Reynolds]

Convert bigWig to WIG via bigWigToWig:

$ bigWigToWig signal.bw signal.wig

Convert WIG to BED via BEDOPS wig2bed:

$ wig2bed < signal.wig > signal.bed

Filter BED by its score column with awk:

$ awk '($5 < 1)' signal.bed > signal.filtered.bed