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 :)
•
written
by
sboardman •
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
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.
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
Hi Jorge, that's great. I wrote a python script to concatenate my covered regions, but this looks much more efficient.
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
Log in to answer this question.
