Bedtools genomecov to identify regions at any coverage
3
0
Entering edit mode
6.8 years ago
sboardman ▴ 10

Hi,

 

Is it possible to use bedtools genomecov with both -bga and -max flags?

I've been trying to do this (code below), but my resulting output isn't binned.

bedtools genomecov -max 1 -bga -ibam input.bam -g hg19.genome > output.bed

chr1  0       554304  0
chr1  554304  554309  5
chr1  554309  554313  6
chr1  554313  554314  1

Whereas I'd like it to be:

chr1  0       554304  0
chr1  554304  554314  1

Any help is much appreciated.

 

SB :)

bedtools genomecov • 7.6k views
ADD COMMENT
2
Entering edit mode
6.8 years ago

The help for the -max option indicates that this does not apply to the bedgraph output.  I would just use awk to "correct" the output:

bedtools genomecov -bga -ibam input.bam -g hg19.genome \
| awk '{if($4>1){$4=1}; print $0}' > output.bed
ADD COMMENT
0
Entering edit mode

Thanks Aaron, that makes sense. Knew I was missing something obvious. I'll work out another way to concatenate the covered regions so I don't have lots of sequential regions with 1 coverage.

ADD REPLY
5
Entering edit mode
6.8 years ago

the -bga option reports depth in BedGraph format including regions with zero coverage. but if you would only need the regions with at least 1 read you can use the -bg option to get them, and ultimately merge them:

bedtools genomecov -ibam input.bam -bg | bedtools merge -i - > covered.bed

you can then get the regions without any reads just by applying a negative search:

awk 'OFS="\t"{print $1,0,$2}' human_hg19.fa.fai | bedtools subtract -a - -b covered.bed > noncovered.bed
ADD COMMENT
0
Entering edit mode

Hi Jorge, that's great. I wrote a python script to concatenate my covered regions, but this looks much more efficient.

ADD REPLY
1
Entering edit mode
6.8 years ago

You could make something like this :

bedtools genomecov -bga -ibam input.bam -g hg19.genome | awk '{if($4<=1){print $0}}' > output.bed

not tested but should work

ADD COMMENT

Login before adding your answer.

Traffic: 1573 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