1
$\begingroup$

I am using bedtools coverage to compute the sequencing depth at every positions of a chromosome but it didn't work as I expected. Instead it reported 0 coverage at every positions.

This is how I did it :

I made a chr1.bed bedfile :

chr start   end 
1   1   158534110

I ran this command :

coverageBed -a -d chr1.bed -b aln.bam -s > coverage.perposition_chr1.tsv

However I don't find any position in the output that is actually covered :

chr start end 1-pos depth
1   1   158534110   1   0
1   1   158534110   2   0
1   1   158534110   3   0
1   1   158534110   4   0
1   1   158534110   5   0
1   1   158534110   6   0
1   1   158534110   7   0
1   1   158534110   8   0
(...)
1   1   158534110   158534109   0

I must do something wrong because I already computed depth at specific pattern positions using the same aln.bam and it was working fine...

Suggestions and tries :

  • Swapping -a & -d parameters produces the same results
  • Adding a strand column to the bedfile still gives 0 coverage.

chr1.bed

1       1       158534110       .       0       +
1       1       158534110       .       0       -

output


1       1       158534110       .       0       +       1       0
1       1       158534110       .       0       +       2       0
(...)
1       1       158534110       .       0       -       158534108       0
1       1       158534110       .       0       -       158534109       0
```
Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ In your .bam file (and the original reference you used) are the chromosomes called '1' '2'; etc, or 'Chr1', 'Chr2' etc? Also, you could try swapping the -a and -d - maybe the -a needs to be immediately before the .bed file. $\endgroup$
    – Jay Moore
    Nov 12, 2019 at 17:18
  • $\begingroup$ @JonathanMoore I bet you're right. The -a expects a string, so that would take -d as the value for -a. Probably worth an answer already. $\endgroup$
    – terdon ♦
    Nov 12, 2019 at 17:21
  • $\begingroup$ @JonathanMoore Yes this is the right chromosome name. I tried swipping parameters and it did produce the same file. @terdon wouldn't I expect an error If it indeed took -d as the -a value ? It is not fitting the right format $\endgroup$
    – Paul Endymion
    Nov 12, 2019 at 17:42
  • 1
    $\begingroup$ Consider posting 1-2 alignments in sam format, which you would expect to give non-zero coverage as the minimum reproducible example. What is the bedtools version? What happens if you drop -s? $\endgroup$
    – Timur Shtatland
    Nov 12, 2019 at 22:52
  • 1
    $\begingroup$ @TimurShtatland Dammit that was it, the bam went blank at some point in my testing. I should have suspected that. It works now. My bad and thanks all for the suggestions.. $\endgroup$
    – Paul Endymion
    Nov 13, 2019 at 9:53

1 Answer 1

Reset to default
2
$\begingroup$

If you prefer to have a more concise report than a per-base one, bedtools genomecov could be a better choice.

Its -bga option allows you to report consecutive bases that have the same coverage as a range in a single line.

In your case, you could achieve this (unfortunately they only offer the strandedness option for a specified strand) by running the following command:

bedtools genomecov -bga -ibam aln.bam -g chr1.bed -strand + > coverage.perposition_chr1.tsv

Output will look like:

chr1  554304  554309  5
chr1  554309  554313  6
chr1  554313  554314  1
chr1  554315  554316  6
chr1  554316  554317  5
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.