Bedtools count split reads spanning whole feature

Question

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

Answer 1

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 2

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