compute coverage with file.bam and bed : bedtools empty output

Question

[noeD]

Hello! :)

I would like to compute the coverage for my file.bam using the information of the file.bed.

My file bed is like this:

chr1    12594   12721   uc010nxq.1_cds_1_0_chr1_12595_f 0   +
chr1    13402   13639   uc010nxq.1_cds_2_0_chr1_13403_f 0   +
chr1    69090   70008   uc001aal.1_cds_0_0_chr1_69091_f 0   +

And I would like to have something like that:

chr1      12594    12721    40
chr1      13402    13639    68

where the IV column is the coverage.

The read in my bam file (I have converted it to sam format in order to check if the chromosome ids in the bam file do not exactly match in the bed file) are like that:

H-D00:95:C62BPANXX:1:2315:2019:91137    163 chr1    10002   38  4S77M1I43M  =   10114   237 CCCAAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTTAACCCTAACCCTAACCCTACCCCTAACCCAAACCTAACCCTA   A@?BCEGFFGDFGFFGDFGFFGDFGFEGDFGFFGDEGFFG@FGFFFDEGEFF@F/BAG/5/EB4;@GE?G/1.BCGAEBFFC3@@GBBDDFGC04<FGF##########################   MC:Z:125M   BD:Z:MMKNOGPPLMNNNOKLMMNOKLMLMNJKLLMNJKLLMNJKLLMNJKLLMNJKLLMNJKLLMOKKLLMOKLMMNOKLMLMOJLLLMMNJLMMNOKLMMNOMMMMMMMMMMMMMNNNNNNOOPMMMM  MD:Z:97A11C10   -PG:Z:MarkDuplicates    RG:Z:HTG-1  BI:Z:QQMTRKRRMSSOQPLRROPPLRROPPKRQOPPLRROPPLRROPPLRROPPLRROPPLSRPPQMSSPQQMSSPQQMSSPQQMSQQPQQMSSPQQMSTQRRRRRRRRRRRRRSSRSSSSSSTUQQQQ  NM:i:3  MQ:i:38 AS:i:103    XS:i:97

I have read this post Tools To Calculate Average Coverage For A Bam File? and first I have tried to compute it with

genomeCoverageBed

but it gave the results for each base.

After that I have tried the following tool and command:

bedtools coverage -a file.bed -b alignment.bam > output

but nothing happened, the output file was empty.

Do you have any suggestions?

Thank you in advance

[Carlo Yague]

What bedtools version are you using ? The coverage command has changed a lot over the time and there can be big differences between the behaviour described in the manual for the last version and the behaviour of a previous version. You could also read the usage with bedtools coverage -h.

[noeD]

thank you for your reply.

I am using this version: 2.25

Answer 1

[Devon Ryan]

bedtools intersect -c will do what you want.

[noeD]

Thank you for your help. I have tried what you suggested:

bedtools intersect -c -abam alignment.bam -b file.bed -bed

These are the first lines of the output:

chr1    10001   10121   H-D00:95:C62BPANXX:1:2315:2019:91137/2  38  +   10001   10121   0,0,0   1   120,    0,  0
chr1    10030   10149   H-D00:95:C62BPANXX:1:2313:14093:20198/2 40  +   10030   10149   0,0,0   1   119,    0,  0

I didn't understand why I tried in my output regions that are not present in my bed file, furthermore I didn't understand the second part of the output... for example in the first line of the output:

+   10001   10121   0,0,0   1   120,    0,  0

Thank you in advance

[Devon Ryan]

Try making the BAM file -b.

[noeD]

I have tried but the file was again empty :(