Problem using bedtools coverage -hist for plotting

[Frederic Thieme]

Hello,

I have performed a targeted sequencing approach and have been trying to plot the coverage of my samples over my target areas. From preliminary analyses I know that my mean coverage is >200 with >90% of the bases (~50 kb in total) covered.

I am using bedtools v2.25.0 on an ubuntu machine. An update to the most current version is not possible at this time.

To do the plotting, I used the protocol given here:

https://github.com/arq5x/bedtools-protocols/blob/master/bedtools.md#ap2b-assessing-coverage-in-exome-capture-experiments

However, my problem becomes apparent at this step:

gcov = cov[cov[,1] == 'all',]

Thus, the output of my analysis only contains two lines for 'all' (0 and 1), suggesting that no coverage greater than 1x was achieved.

To make sure that this was not due to a fault in my data I decided to follow the tutorial using the exome .bam file from 1000 genomes (also found in the BEDtools paper from 2015, basic protocol 3).

I used the following command to perform the analysis:

bedtools coverage -hist -abam NA12891.exome.bam -b targets.numeric.chroms.bed > NA12891.exome.coverage.hist.txt

Here are the first three lines from the resulting file (the last column does not look as it does in the paper):

1       10004   10078   SRR098359.28989051/1    0       +       10004   10078   0,0,0   1       74,0,       0       74      74      1.0000000
1       10005   10043   SRR098359.37775242/2    0       +       10005   10043   0,0,0   1       38,0,       0       38      38      1.0000000
1       10014   10090   SRR098359.113245777/2   0       -       10014   10090   0,0,0   1       76,0,       0       76      76      1.0000000

And the last three lines:

hs37d5  35474151        35474197        SRR098359.111788980/2   13      +       35474151        35474197    0,0,0   1       46,     0,      0       46      46      1.0000000
all     0       794493821       1572993857      0.6490867
all     1       778500036       1572993857      0.3509133

Again, only two different coverages (0x and 1x) are given, while I expected several more (also according to the paper).

I did not plot this; plotting the result of my previous attempt using my own data showed a plot with a straight line between the two generated points (0x and 1x coverage, respectively).

What am I doing wrong?

Thank you and best wishes,

Frederic

[Aaron Quinlan]

Hi Frederic,

The problem you are facing is caused by the fact that the usage of coverage has changed.  

See “Important” here: http://bedtools.readthedocs.io/en/latest/content/tools/coverage.html

For your use case, you want the BED file to be -a and the BAM file to be -b and you want both files to be sorted by chromosome and then by start coordinate.  For example,

sort -k1,1 -k2,2n targets.numeric.chroms.bed > targets.numeric.chroms.sorted.bed

samtools sort NA12891.exome.bam > NA12891.exome.sorted.bam\

bedtools coverage -hist -a targets.numeric.chroms.sorted.bed -b NA12891.exome.sorted.bam > NA12891.exome.coverage.hist.txt

Also, I would recommend updating to version 2.27.1, as bugs and performance improvements have been applied to this tool.

Aaron

--
You received this message because you are subscribed to the Google Groups "bedtools-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to bedtools-discu...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.