Bedtools Coverage 1x Max Issue

Question

[kerriganblake]

I am trying to make a histogram of the sequencing depth in a set of BAM files using the WES kit's target_regions.bed file. I ran the following command with no error messages generated:

bedtools coverage -hist -abam sample.bam -b target_regions.bed | grep ^all > sample.bam.hist.all.txt

However, after running this on each of my 35 samples, all of the files show at most 1x coverage (which can be clearly demonstrated as wrong in IGV), like the following:

all     0       1644941785      -277541545      0.4094517
all     1       -1922483330     -277541545      0.5905483

I did ensure the BED file was sorted, though I noticed that the position of chrM is in a different location in my BAM and my BED. Could this cause the issue or is it something more complicated? I've been reading forums all day to no avail, so any suggestions are welcome! Thanks!

[stachele]

Did you make any progress on this problem? Currently googling and having trouble making any headway...

Answer 1

[melissawongukm]

Hi, it seems like the bedtools command has changed. The correct command should be:

bedtools coverage -hist -a x.bed -b x.bam

Now the file should contain each feature on the bed file plus "all" lines. To make sure all lines contain 7 columns and can be read by R read.table, I added two more columns to "all" lines.

bedtools coverage -hist -a x.bed -b x.bam | awk '{if ($1=="all") print $0"\t.\t."; else print $0}' > x.hist.txt