I have a BAM file produced using BWA-MEM and I need to calculate the coverage of each exon in this sequencing. I tried using this command:
bedtools coverage -a exome_sorted.bed -b exome_alignment.bam -sorted
And the result was the following:
chr1 12595 12802 ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242 . ref|DDX11L1 0 0 207 0.0000000 chr1 13163 13658 ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242 . ref|DDX11L1 0 0 495 0.0000000 chr1 14620 15015 ref|WASH7P,ref|NR_024540,ens|ENST00000538476,ens|ENST00000438504,ens|ENST00000488147,ens|ENST00000541675,ens|ENST00000423562 . ref|WASH7P 0 0 395 0.0000000 ...
That is, zero coverage for all exons. But there are reads mapping to the exons, as I can see by opening my BAM file on IGV and looking at some exonic regions:
I have already tried sorting the BED file using
sort -k1,1 -k2,2n (as recommended on the Bedtools documentation) and using
bedtools sort with the
-faidx option, but neither option solved my issue.
if I skip the
-sorted option, the memory usage becomes way too high.
I have also checked if the chromosome names are the same in the BAM file and in the BED file, using
samtools view -H. Both use
chr1,chr2,...,chr22,chrX,chrY, although BAM has
chrM and stuff like
chrUn_gl000214 which are not on the BED.
@ATpoint pointed out that all visible reads in this region have MAPQ0. I looked at another exon that did not have this issue, but Bedtools still outputs zero coverage. IGV picture below: