Bedtools genomecov to identify regions at any coverage

Question

[sboardman]

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 :)

Answer 1

[Aaronquinlan]

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

[sboardman]

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.

Answer 2

[Jorge Amigo]

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

[sboardman]

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

Answer 3

[Nicolas Rosewick]

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