3
$\begingroup$

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:

IGV view of the exonic region chr1:12,595-12,802 , clearly showing many reads mapping.

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.

Thanks!

EDIT:

@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:

IGV view of the exonic region chr1:12,595-12,802 , clearly showing many reads mapping.

$\endgroup$
0
2
$\begingroup$

If this is your output :

chr1    12595   12802   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       207     0.0000000

Then I guess your bedfile looks like this (the first 6 fields) :

  • chrom : chr1
  • FeatureStart : 12595
  • FeatureEnd : 12802
  • name : ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242
  • score : .
  • strand : ref|DDX11L1

The issue I see here is that you are giving the wrong information to the strand field of your bedfile. It should be either +, either - or . if you don't want to consider the strand information.

Hope it helps !

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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