text: user9497: [Note: this is for bedtools v2.30]
  I want to count the number of RNA-Seq reads that fully cover a given splice junction. For that, I thought of defining a BED feature around the junction, then using bedtools coverage.
  Here is a simplified situation:
  1-BASED POS 1|      10|       20|       30|
  CHROMOSOME   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  BED FILE              **                   
  BAM READ FWD >>>>>>>>>>                    
  BAM READ REV                  <<<          
  
  We can generate the corresponding example files:
  echo -e "I\t9\t11\tmy_position\t.\t+" > single_place.bed
  
  echo -e "@HD\tVN:1.4\tSO:coordinate
  @SQ\tSN:I\tLN:30
  first_read\t99\tI\t1\t255\t10M\t=\t20\t20\tAATTCCGGAA\tDDDDDDDDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0
  first_mate\t147\tI\t20\t255\t3M\t=\t1\t20\tATC\tDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0" | samtools view -b > fwd_only.bam
  samtools index fwd_only.bam
  
  
  cat single_place.bed
  I       9       11      my_position     .       +
  samtools view fwd_only.bam
  first_read      99      I       1       255     10M     =       20      20      AATTCCGGAA      DDDDDDDDDD      NH:i:1  HI:i:1  AS:i:127        nM:i:0
  first_mate      147     I       20      255     3M      =       1       20      ATC     DDD     NH:i:1  HI:i:1  AS:i:127        nM:i:0
  
  The per-base count (with the -d option) gives the expected result:
  coverageBed -d -a single_place.bed -b fwd_only.bam
  I       9       11      my_position     .       +       1       1
  I       9       11      my_position     .       +       2       0
  
  And using a filter -f 1.0 to only keep reads that cover the entire feature indeed removes our read:
  coverageBed -a single_place.bed -b fwd_only.bam
  I       9       11      my_position     .       +       1       1       2       0.5000000
  coverageBed -f 1.0 -a single_place.bed -b fwd_only.bam
  I       9       11      my_position     .       +       0       0       2       0.0000000
  
  However, in practice, I have a lot of spliced reads, so I wish to use the -split option. In that case, the filtering doesn't work anymore:
  coverageBed -split -f 1.0 -a single_place.bed -b fwd_only.bam
  I       9       11      my_position     .       +       1       1       2       0.5000000
  
  even though the per-base count is still unchanged:
  coverageBed -d -split -a single_place.bed -b fwd_only.bam
  I       9       11      my_position     .       +       1       1
  I       9       11      my_position     .       +       2       0
  
  Why does -split interfere with -f 1.0?
  With spliced reads
  Note, in practice I'm more interested in cases where the reads can be spliced, such as this:
  1-BASED POS 1|      10|       20|       30|
  CHROMOSOME   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  BED FILE              **                    
  BAM READ FWD >>>>>>>>>>__>>                
  BAM READ REV                  <<<          
  
  where _ is a spliced read. The corresponding BAM can be generated with:
  echo -e "@HD\tVN:1.4\tSO:coordinate
  @SQ\tSN:I\tLN:30
  first_read\t99\tI\t1\t255\t10M2N2M\t=\t20\t20\tAAAAAAAAAACC\tDDDDDDDDDDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0
  first_mate\t147\tI\t20\t255\t3M\t=\t1\t20\tATC\tDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0" | samtools view -b > fwd_only.bam
  samtools index fwd_only.bam
  
  In that situation, it seems I need to use -split and -f 1.0. I thought of filtering a posteriori with e.g. awk '$10>0.5, but that doesn't seem to work either: when not using -d, the output of coverageBed does not separate the reads that cover the feature partially vs fully.
  
  user73: Not quite sure if this is the issue you're having, but I've encountered a similar situation with coverage estimation when trying to map nanopore cDNA reads to the genome. I have a full workflow for mapping oriented long reads here, but here's a summary of the split issue:
  
  The default options for BEDTools treat sequence deletions (which happen frequently in nanopore reads) as a drop in coverage, which can make exon hunting and coverage calculation more difficult. I submitted a pull request to add an additional ignoreD parameter to the command line to allow cDNA reads with split coverage across introns to ignore deletions when considering coverage; this request has now been incorporated into the main BEDtools repository (as of v2.30.0).
  
  Depending on how your mapping is set up, adding -ignoreDto your command line may allow spliced reads to be properly counted for coverage with no deletion gaps:
  $ bedtools genomecov -bga -strand '+' -split -ignoreD -ibam mapped_BC01.sam
  chrX    0       7384244 0
  chrX    7384244 7384413 1
  chrX    7384413 171031299       0
  
  $ bedtools genomecov -bga -strand '+' -split -ignoreD -ibam mapped_BC01_exon.sam
  chr11   0       62551499        0
  chr11   62551499        62551595        1
  chr11   62551595        62552139        0
  chr11   62552139        62553212        1
  chr11   62553212        122082543       0
  
  
  user9497: It appears this is an issue that has been raised on Github, and that I had unfortunately missed.
  This issue happens in bedtools coverage (my examples above), but also in bedtools intersect. So, I believe this is actually the same issue as #928. Indeed, switching to bedtools v2.27 solves my issue described above (with v2.30).
  Note that this particular issue likely appeared as a consequence of solving issues #773, #750, #673.
  My workaround is to use bedtools intersect INVERTING the bed and bam, so that I just filter in the reads that match the feature, then count the resulting number of rows.
  intersectBed -split -bed -wb -F 1.0 \
               -a fwd_only.bam -b single_place.bed | \
    cut -f16 | \
    sort | \
    uniq -c
  
  
original_text: 9497: <p>[Note: this is for bedtools v2.30]</p>
  <p>I want to count the number of RNA-Seq reads that fully cover a given splice junction. For that, I thought of defining a BED feature around the junction, then using <code>bedtools coverage</code>.</p>
  <p>Here is a simplified situation:</p>
  <pre><code>1-BASED POS 1|      10|       20|       30|
  CHROMOSOME   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  BED FILE              **                   
  BAM READ FWD &gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;                    
  BAM READ REV                  &lt;&lt;&lt;          
  </code></pre>
  <p>We can generate the corresponding example files:</p>
  <pre><code>echo -e &quot;I\t9\t11\tmy_position\t.\t+&quot; &gt; single_place.bed
  
  echo -e &quot;@HD\tVN:1.4\tSO:coordinate
  @SQ\tSN:I\tLN:30
  first_read\t99\tI\t1\t255\t10M\t=\t20\t20\tAATTCCGGAA\tDDDDDDDDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0
  first_mate\t147\tI\t20\t255\t3M\t=\t1\t20\tATC\tDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0&quot; | samtools view -b &gt; fwd_only.bam
  samtools index fwd_only.bam
  
  
  cat single_place.bed
  I       9       11      my_position     .       +
  samtools view fwd_only.bam
  first_read      99      I       1       255     10M     =       20      20      AATTCCGGAA      DDDDDDDDDD      NH:i:1  HI:i:1  AS:i:127        nM:i:0
  first_mate      147     I       20      255     3M      =       1       20      ATC     DDD     NH:i:1  HI:i:1  AS:i:127        nM:i:0
  </code></pre>
  <p>The per-base count (with the <code>-d</code> option) gives the expected result:</p>
  <pre><code>coverageBed -d -a single_place.bed -b fwd_only.bam
  I       9       11      my_position     .       +       1       1
  I       9       11      my_position     .       +       2       0
  </code></pre>
  <p>And using a filter <code>-f 1.0</code> to only keep reads that cover the entire feature indeed removes our read:</p>
  <pre><code>coverageBed -a single_place.bed -b fwd_only.bam
  I       9       11      my_position     .       +       1       1       2       0.5000000
  coverageBed -f 1.0 -a single_place.bed -b fwd_only.bam
  I       9       11      my_position     .       +       0       0       2       0.0000000
  </code></pre>
  <p>However, in practice, I have a lot of spliced reads, so I wish to use the <code>-split</code> option. In that case, the filtering doesn't work anymore:</p>
  <pre><code>coverageBed -split -f 1.0 -a single_place.bed -b fwd_only.bam
  I       9       11      my_position     .       +       1       1       2       0.5000000
  </code></pre>
  <p>even though the per-base count is still unchanged:</p>
  <pre><code>coverageBed -d -split -a single_place.bed -b fwd_only.bam
  I       9       11      my_position     .       +       1       1
  I       9       11      my_position     .       +       2       0
  </code></pre>
  <p>Why does <code>-split</code> interfere with <code>-f 1.0</code>?</p>
  <h3>With spliced reads</h3>
  <p>Note, in practice I'm more interested in cases where the reads can be spliced, such as this:</p>
  <pre><code>1-BASED POS 1|      10|       20|       30|
  CHROMOSOME   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  BED FILE              **                    
  BAM READ FWD &gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;__&gt;&gt;                
  BAM READ REV                  &lt;&lt;&lt;          
  </code></pre>
  <p>where <code>_</code> is a spliced read. The corresponding BAM can be generated with:</p>
  <pre><code>echo -e &quot;@HD\tVN:1.4\tSO:coordinate
  @SQ\tSN:I\tLN:30
  first_read\t99\tI\t1\t255\t10M2N2M\t=\t20\t20\tAAAAAAAAAACC\tDDDDDDDDDDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0
  first_mate\t147\tI\t20\t255\t3M\t=\t1\t20\tATC\tDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0&quot; | samtools view -b &gt; fwd_only.bam
  samtools index fwd_only.bam
  </code></pre>
  <p>In that situation, it seems I need to use <code>-split</code> and <code>-f 1.0</code>. I thought of filtering <em>a posteriori</em> with e.g. <code>awk '$10&gt;0.5</code>, but that doesn't seem to work either: when not using <code>-d</code>, the output of <code>coverageBed</code> does not separate the reads that cover the feature partially vs fully.</p>
   73: <p>Not quite sure if this is the issue you're having, but I've encountered a similar situation with coverage estimation when trying to map nanopore cDNA reads to the genome. I have a full workflow for mapping oriented long reads <a href="https://dx.doi.org/10.17504/protocols.io.b54iq8ue" rel="nofollow noreferrer">here</a>, but here's a summary of the split issue:</p>
  <blockquote>
  <p>The default options for BEDTools treat sequence deletions (which happen frequently in nanopore reads) as a drop in coverage, which can make exon hunting and coverage calculation more difficult. I submitted <a href="https://github.com/arq5x/bedtools2/pull/868" rel="nofollow noreferrer">a pull request</a> to add an additional ignoreD parameter to the command line to allow cDNA reads with split coverage across introns to ignore deletions when considering coverage; this request has now been incorporated into the main BEDtools repository (as of v2.30.0).</p>
  </blockquote>
  <p>Depending on how your mapping is set up, adding <code>-ignoreD</code>to your command line may allow spliced reads to be properly counted for coverage with no deletion gaps:</p>
  <pre><code>$ bedtools genomecov -bga -strand '+' -split -ignoreD -ibam mapped_BC01.sam
  chrX    0       7384244 0
  chrX    7384244 7384413 1
  chrX    7384413 171031299       0
  
  $ bedtools genomecov -bga -strand '+' -split -ignoreD -ibam mapped_BC01_exon.sam
  chr11   0       62551499        0
  chr11   62551499        62551595        1
  chr11   62551595        62552139        0
  chr11   62552139        62553212        1
  chr11   62553212        122082543       0
  </code></pre>
   9497: <p>It appears this is an issue that has been raised on Github, and that I had unfortunately missed.</p>
  <p>This issue happens in <code>bedtools coverage</code> (my examples above), but also in <code>bedtools intersect</code>. So, I believe this is actually the same issue as <a href="https://github.com/arq5x/bedtools2/issues/928" rel="nofollow noreferrer">#928</a>. Indeed, switching to bedtools v2.27 solves my issue described above (with v2.30).</p>
  <p>Note that this particular issue likely appeared as a consequence of solving issues #773, #750, #673.</p>
  <p>My workaround is to use <code>bedtools intersect</code> INVERTING the bed and bam, so that I just filter in the reads that match the feature, then count the resulting number of rows.</p>
  <pre><code>intersectBed -split -bed -wb -F 1.0 \
               -a fwd_only.bam -b single_place.bed | \
    cut -f16 | \
    sort | \
    uniq -c
  </code></pre>
  