2
$\begingroup$

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

$\endgroup$

2 Answers 2

1
$\begingroup$

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
$\endgroup$
1
$\begingroup$

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
$\endgroup$
1
  • 1
    $\begingroup$ I do not think this solves it, as the split is to handle the Ns 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. $\endgroup$
    – Alexlok
    Mar 9 at 23:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.