Bedtools coverage gives zero coverage output

Question

[biotech]

I'm running bedtools coverage to get coverage for each gene but I have zero coverage output in resulting bed file. I think the problem is in the bed file since is a custom one.

Here are my files and commands:

bedtools coverage -abam MZ123.bam -b 7a1.bed > coverage.txt

coverage.txt

cpdna   1   865 Ga0114182_111   1000    +   1   865 65,105,225  0   0   864 0.0000000
cpdna   841 2053    Ga0114182_112   1000    -   841 2053    65,105,225  0   0   1212    0.0000000
cpdna   2174    2759    Ga0114182_113   1000    -   2174    2759    65,105,225  0   0   585 0.0000000
cpdna   2755    3886    Ga0114182_114   1000    -   2755    3886    65,105,225  0   0   1131    0.0000000
cpdna   4008    4587    Ga0114182_115   1000    -   4008    4587    65,105,225  0   0   579 0.0000000
cpdna   5047    5419    Ga0114182_116   1000    +   5047    5419    65,105,225  0   0   372 0.0000000
cpdna   5681    6212    Ga0114182_117   1000    +   5681    6212    65,105,225  0   0   531 0.0000000
cpdna   6324    7119    Ga0114182_118   1000    -   6324    7119    65,105,225  0   0   795

And how my input bed looks:

track name=vitVinGenes description="V. vinifera cpdna genes" itemRgb=On
cpdna   1   865 Ga0114182_111   1000    +   1   865 65,105,225
cpdna   841 2053    Ga0114182_112   1000    -   841 2053    65,105,225
cpdna   2174    2759    Ga0114182_113   1000    -   2174    2759    65,105,225
cpdna   2755    3886    Ga0114182_114   1000    -   2755    3886    65,105,225
cpdna   4008    4587    Ga0114182_115   1000    -   4008    4587    65,105,225
cpdna   5047    5419    Ga0114182_116   1000    +   5047    5419    65,105,225
cpdna   5681    6212    Ga0114182_117   1000    +   5681    6212    65,105,225
cpdna   6324    7119    Ga0114182_118   1000    -   6324    7119    65,105,225
cpdna   7153    7771    Ga0114182_119   1000    -   7153    7771    65,105,225
cpdna   7767    8382    Ga0114182_1110  1000    -   7767    8382    65,105,225
cpdna   8381    9401    Ga0114182_1111  1000    -   8381    9401    65,105,225
cpdna   9537    11523   Ga0114182_1112  1000    +   9537    11523   65,105,225
cpdna   11632   12187   Ga0114182_1113  1000    +   11632   12187   65,105,225
cpdna   12304   15025   Ga0114182_1114  1000    +   12304   15025   65,105,225

Some more hints. According to flagstat, reads were mapped:

samtools flagstat MZ123.bam
        1367199 + 0 in total (QC-passed reads + QC-failed reads)
        1259488 + 0 mapped (92.12% : N/A)

[James Ashmore]

Your bam file should be the -b file and your bed file should be the -a file. The coverage command reports the number of features in B that overlap the the A interval. In your case that would be the number of aligned reads in the BAM file which overlap each gene in the BED file.

[biotech]

This command you refer is creating a HUGE output. Is this what is expected? Seems command is not finishing. Something is wrong here. I just have about 5K features.

wc -l coverage.txt
72317662

and still running...

[James Ashmore]

Which version of bedtools are you using?

[biotech]

I'm using bedtools v2.24.0 for mac

[James Ashmore]

According to the release history, v2.24.0 coverage command should take the intervals as the a file, and the alignments as the b file. Try updating to the latest version. Are the chromosome names you are using in the BED file the same as those found in the BAM file?

[biotech]

No, this is wrong. I will fix chromosome names and then run the way you say.