Question: Bedtools genomecov to identify regions at any coverage
0
gravatar for sboardman
4.4 years ago by
sboardman0
United Kingdom
sboardman0 wrote:

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

genomecov bedtools • 5.0k views
ADD COMMENTlink modified 4.4 years ago by Jorge Amigo11k • written 4.4 years ago by sboardman0
2
gravatar for Aaronquinlan
4.4 years ago by
Aaronquinlan11k
United States
Aaronquinlan11k wrote:

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 COMMENTlink written 4.4 years ago by Aaronquinlan11k

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 REPLYlink written 4.4 years ago by sboardman0
5
gravatar for Jorge Amigo
4.4 years ago by
Jorge Amigo11k
Santiago de Compostela, Spain
Jorge Amigo11k wrote:

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 COMMENTlink modified 4.4 years ago • written 4.4 years ago by Jorge Amigo11k

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

ADD REPLYlink written 4.4 years ago by sboardman0
0
gravatar for Nicolas Rosewick
4.4 years ago by
Belgium, Brussels
Nicolas Rosewick7.6k wrote:

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 COMMENTlink written 4.4 years ago by Nicolas Rosewick7.6k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1703 users visited in the last hour