What is the correct way to get the total number of reads strictly contained in each interval in a GFF from a BAM file while enforcing strandedness? What I am looking for is very close to this intersectBed feature:
-c For each entry in A, report the number of overlaps with B.
- Reports 0 for A entries that have no overlap with B.
- Overlaps restricted by -f and -r.
Except that I'd like the number of overlaps in A for each entry in B (i.e. the other way around). If I do:
intersectBed -abam mybam.bam -b mygff.gff -s -f 1 -wb
Then my understanding is that this will report the entry in B for each overlap with A. But I'd like each entry in B to be outputted exactly once, with the number of reads from A that are contained strictly within it. I'm not sure how to enforce strict containment here.
Is coverageBed the solution to this? Or multicov? I'm not sure how to enforce strict containment using coverageBed - it's not clear to me if that's the default from the docs. Thanks.
•
written
by
user
2 answers
Another way to do this is to use the bedmap application from the BEDOPS suite for this, where your intervals are reference regions and reads are the map regions.
First, convert the GFF to sorted BED with the gff2bed script and sort-bed:
$ gff2bed.py < mygff.gff | sort-bed - > myGff.bed
Convert the reads from BAM to sorted BED with the bam2bed script:
$ bam2bed.csh < mybam.bam | sort-bed - > myBam.bed
With bedmap, add the --echo operator to output each reference interval from myGff.bed. Add the --fraction-map 1 operator to enforce the overlap criteria that mapped reads must be contained entirely within a reference interval. Add the --count operator to yield the number of mapped reads within each reference interval, which are printed adjacent to the reference interval data.
In sum, as an example:
$ bedmap --echo --fraction-map 1 --count myGff.bed myBam.bed > myAnswer.bed
If you need to apply strandedness criteria, such as shifting reverse-strand reads by a base, the bedops --range L:R operator can be used to pre-adjust element coordinates asymmetrically, e.g. bedops --range -1:-1 --everything foo.bed > bar.bed will reduce each start and stop coordinate by one base. Here, you might split your inputs by strand, pre-adjust a stranded input, run two bedmap operations on each input and then collate the two results at the end.
•
written
by
Alex Reynolds
No problem. Hopefully it helps someone else trying to do the same thing.
Currently, the coverage tool cannot enforce fractional coverage. However, one strategy would be to use intersect to enforce the strict containment and then pass only those alignments that meet this criteria to the coverage tool. The example below demonstrates this strategy.
bedtools intersect -abam mybam.bam -b mygff.gff -s -f 1 -ubam | \
bedtools coverage -abam - -b mygff.gff
Perhaps look into mkfifo and named pipes: http://www.linuxjournal.com/article/2156
Log in to answer this question.
-
bedtools intersectBed and GRanges countOverlaps
written by User 7754
-
Question about number of reads within intervals
written by iraun
-
How To Count Stand-Specific Paired-End Rna-Seq Reads Overlapping Known Protein Coding Genes ?
written by biorepine
-
Bowtie2 Mapping Different Number Of Reads To Same Sequence When Ref-Seq Is Part Of Different Indexes
written by rohita.sinha •
-
Getting All Reads That Align To A Region In Compact Bed Format Using Bedtools?
written by user
-
Intersectbed/Coveragebed -Split Purify Exon?
written by Puriney
-
How To Get All Entries Of B With Bedtools
written by lyfsa •
-
Reporting The Bam Reads Overlapping A Set Of Intervals With Bedtools
written by User 9996
-
Bedtools Compare Multiple Bed Files?
written by Bioscientist
-
Counting Number Of Bam Reads Directly Within Set Of Intervals With Bedtools
written by User 9996
