Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-03-10T00:16:43.427Z
Answers.LastActivityDate: 2022-03-10T00:16:43.427Z
Answers.UserDisplayName: Alexlok
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-03-09T23:51:31.460Z
Answers.Text: I do not think this solves it, as the `split` is to handle the `N`s in the CIGAR string, while these are short reads with no deletion. But by pointing to your PR, you made me find a relevant Github issue that I had somehow missed, thank you! I'll post it as an answer for future reference.
Answers.UserDisplayName: Alexlok
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-03-08T22:45:25.240Z
Answers.LastActivityDate: 2022-03-10T03:03:37.427Z
Answers.LastEditDate: 2022-03-10T03:03:37.427Z
Answers.UserDisplayName: gringer
Body: <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>
  
Tags: <bedtools>
Title: Bedtools count split reads spanning whole feature
ThreadText: # Bedtools count split reads spanning whole feature
  
  (Asked by: username_0 on 2022-03-07)
  
  [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.
  
  **Tags:** bedtools
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_0 on 2022-03-10. Score: 1 | Accepted answer
  
  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](https://github.com/arq5x/bedtools2/issues/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
  
  ```
  
  
  ---
  
  ### Answer by username_1 on 2022-03-08. Score: 1
  
  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](https://dx.doi.org/10.17504/protocols.io.b54iq8ue), 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](https://github.com/arq5x/bedtools2/pull/868) 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 `-ignoreD`to 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
  
  ```
  
  > **Comments:**
  > **username_0**: I do not think this solves it, as the `split` is to handle the `N`s in the CIGAR string, while these are short reads with no deletion. But by pointing to your PR, you made me find a relevant Github issue that I had somehow missed, thank you! I'll post it as an answer for future reference.  
  
  ---
  
  