Question

Bedtools intersect PE BAM reads

1

Entering edit mode

3.6 years ago

maxrwjones ▴ 60

Hi all, hoping someone can help with a nit-picky question!

When bedtools intersect is provided with a BAM file generated from paired-end reads as the query file (-a option), does it look for an intersection of the individual reads with the reference sequences or does it use the inserts?

e.g. for these paired entries in my BAM (converted to SAM, raw seq omitted) file, would it look for intersections with the feature [1,298,531 to 1,298,889] (i.e. the insert) or with [1,298,531 to 1,298,681] and [1,298,739 to 1,298,889] (i.e. the individual mates). The difference being whether annotations falling within the 'unknown' insert region will count as having been intersected or not.

A00877:614:HWGW7DSXY:1:2331:23204:15123 163     chr1A   1298531 42      150M    =       1298739 358
A00877:614:HWGW7DSXY:1:2331:23204:15123 83      chr1A   1298739 42      150M    =       1298531 -358

Hope I'm making sense!

intersect BAM bedtools • 2.3k views

ADD COMMENT • link 3.6 years ago by maxrwjones ▴ 60

score 1 · Answer 1 · 2021-05-20

From the Bedtools manual:

Note that only the mate in the BAM alignment is compared to the BED file. Thus, if only one end of a paired-end sequence overlaps with a feature in B, then that end will be written to the BAM output. By contrast, the other mate for the pair will not be written. One should use pairToBed(Section 5.2) if one wants each BAM alignment for a pair to be written to BAM output.

going on that logic I would assume it threads (by default) both reads of a pair as separate entities. To me, this sounds that if the region you are querying is in the 'insert' it will not be reported. At first sight you can overcome this with the use of pairToBed or such? (Extract Only Paired-End Reads That Map A Specific Interval)