This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How To Compute The Number Of Reads In Bam That Overlap With A Region? The Overlap Must Be 100%.

how to compute the number of reads in BAM that overlap with a region? the overlap must be 100%.

reads overlap bam

3 answers

It is unclear from your question whether you want a read to be contained entirely within a region, or a region contained entirely within a read, when counting reads.

But that's okay! You can use the --count operator in BEDOPS bedmap with BEDOPS bam2bed to efficiently calculate either result, specifying different overlap criteria in either case.

First, use BEDOPS sort-bed to make sure your input regions are ready for use with bedmap:

$ sort-bed unsortedRegions.bed > sortedRegions.bed

Preparing input with sort-bed allows BEDOPS tools to work fast and have a very low memory profile, compared with alternatives.

We're now ready to search for reads that overlap regions.

Let's examine the first overlap scenario: Use the --fraction-map operator with bedmap to enforce the overlap criterion that a BAM read (a "mapped" element) must overlap the region by 100%, i.e. that a read is contained entirely within a region-of-interest:

$ bam2bed < reads.bam \
    | bedmap --echo --count --fraction-map 1.0 sortedRegions.bed - \
    > answer.bed

The format of answer.bed will be:

[ region-1 ] | [ number of BAM reads entirely within region-1 ]
[ region-2 ] | [ number of BAM reads entirely within region-2 ]
...
[ region-N ] | [ number of BAM reads entirely within region-N ]

Now let's discuss the second overlap case. If, instead, you want the region-of-interest (a "reference" element) to be contained entirely within the read, then use the --fraction-ref operator:

$ bam2bed < reads.bam \
    | bedmap --echo --count --fraction-ref 1.0 sortedRegions.bed - \
    > answer.bed

The format of answer.bed will be:

[ region-1 ] | [ # BAM reads that region-1 is entirely contained within ]
[ region-2 ] | [ # BAM reads that region-2 is entirely contained within ]
...
[ region-N ] | [ # BAM reads that region-N is entirely contained within ]

You can give any fraction from 0 to 1, inclusive, to --fraction-map and --fraction-ref, depending on the stringency of overlap that you need.

Other overlap criteria are available, which can be useful, depending on your analysis.

You can use bedtools :

something like this (maybe someone can correct me, I'm not sure of the params):

intersectBed -abam in.bam -b region.gtf -c -f 1

I think this can only find B overlap A with 100%,its A 's 100%, not B'100%

Do you want to find reads that entirely contain a small region, or reads that are entirely within a region larger than the read?

You can also try htseq

Log in to answer this question.