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

CC BY-SA 4.0
2022-03-10T00:16:43.427Z
2022-03-10T00:16:43.427Z
Alexlok
<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>

CC BY-SA 4.0
2022-03-09T23:51:31.460Z
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.
Alexlok
CC BY-SA 4.0
2022-03-08T22:45:25.240Z
2022-03-10T03:03:37.427Z
2022-03-10T03:03:37.427Z
gringer
<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>

<bedtools>
Bedtools count split reads spanning whole feature
# 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. 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, 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

> **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.  

---

