Biostar Test Site

This is site is used for testing only. Visit: https://www.biostars.org to ask a question.

bigWig to bed for regions above/below threshold
4
1
Entering edit mode
3.6 years ago
skanterakis ▴ 110

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!

bed bigwig • 3.2k views
ADD COMMENT
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

There is no fixed meaning to that value.

ADD REPLY
1
Entering edit mode
3.6 years ago

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.

ADD COMMENT
1
Entering edit mode
3.6 years ago

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

ADD COMMENT
1
Entering edit mode
3.6 years ago
  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.
ADD COMMENT
0
Entering edit mode
3.6 years ago

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
ADD COMMENT

Login before adding your answer.

Traffic: 181 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6