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!
4 answers
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.
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
- Use
import()from the rtracklayer bioconductor package to read the bigwig file. - Use
slice()from the GenomicRanges bioconductor package on the results from #1. - Use
export()from the rtracklayer bioconductor package on the results from #2 to write the BED file.
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
Log in to answer this question.
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.
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
There is no fixed meaning to that value.