text: user3389: 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.
  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:
  
  
  user6273: 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 !
  
original_text: 3389: <p>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:</p>
  <p><code>bedtools coverage -a exome_sorted.bed -b exome_alignment.bam -sorted</code></p>
  <p>And the result was the following:</p>
  <pre><code>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
  ...
  </code></pre>
  <p>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:</p>
  <p><a href="https://i.stack.imgur.com/TdnLY.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/TdnLY.png" alt="IGV view of the exonic region chr1:12,595-12,802 , clearly showing many reads mapping." /></a></p>
  <p>I have already tried sorting the BED file using <code>sort -k1,1 -k2,2n</code> (<a href="https://bedtools.readthedocs.io/en/latest/content/tools/coverage.html" rel="nofollow noreferrer">as recommended on the Bedtools documentation</a>) and using <code>bedtools sort</code> with the <code>-faidx</code> option, but neither option solved my issue.</p>
  <p>if I skip the <code>-sorted</code> option, the memory usage becomes way too high.</p>
  <p>I have also checked if the chromosome names are the same in the BAM file and in the BED file, using <code>samtools view -H</code>. Both use <code>chr1,chr2,...,chr22,chrX,chrY</code>, although BAM has <code>chrM</code> and stuff like <code>chrUn_gl000214</code> which are not on the BED.</p>
  <p>Thanks!</p>
  <p>EDIT:</p>
  <p>@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:</p>
  <p><a href="https://i.stack.imgur.com/qJvWz.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/qJvWz.png" alt="IGV view of the exonic region chr1:12,595-12,802 , clearly showing many reads mapping." /></a></p>
   6273: <p>If this is your output :</p>
  
  <pre><code>chr1    12595   12802   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       207     0.0000000
  </code></pre>
  
  <p>Then I guess your bedfile looks like this (the first 6 fields) :</p>
  
  <ul>
  <li>chrom : chr1</li>
  <li>FeatureStart : 12595</li>
  <li>FeatureEnd : 12802</li>
  <li>name : ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242</li>
  <li>score : .</li>
  <li>strand : ref|DDX11L1</li>
  </ul>
  
  <p>The issue I see here is that you are giving the wrong information to the strand field of your bedfile. It should be either <code>+</code>, either <code>-</code> or <code>.</code> if you don't want to consider the strand information.</p>
  
  <p>Hope it helps !</p>
  