How to extract reads that map within a defined interval from each other

Question

[dgimode]

I have paired end sequence data which I have mapped to a reference. I need to extract reads from the bam file that map either next to each other or within a defined interval from each other, for example, between 200 and 300 bp from each other. Is it possible to do this? For instance, if my first read pair map between position 90,000 and 90,130 and a second read pair maps between 90,380 and 90,500. These 2 reads map 250 bp apart on the reference. Is it possible to extract all the reads that have this characteristic? Thanks

[venu]

Did you check bedtools closestBed

Answer 1

[Alex Reynolds]

Sure, you can do this with a few genomic set operations with BEDOPS tools:

1. Convert your reads to sorted BED with bam2bed
2. Use bedmap --range to find elements within the first, smaller range boundary (set A)
3. Use bedmap --range again to find elements within the second, larger range boundary (set B)
4. Use bedops --not-element-of to filter out elements from set B that overlap elements in set A

What's left at the end are your elements between the ranges specified for sets A and B.

For instance, to get reads that are between 200 and 300 bases from one another:

$ bam2bed < reads.bam > reads.bed
$ bedmap --range 200 --count --echo-map --delim '\t' reads.bed | awk '$1>1' | cut -f2- | tr ';' '\n' | sort-bed - > reads_within_200_bases.bed
$ bedmap --range 300 --count --echo-map --delim '\t' reads.bed | awk '$1>1' | cut -f2- | tr ';' '\n' | sort-bed - > reads_within_300_bases.bed
$ bedops --not-element-of 1 reads_within_300_bases.bed reads_within_200_bases.bed > answer.bed

The file answer.bed has reads between 200 and 300 bases of one another.

Answer 2

[dgimode]

The BEDOPS way worked just fine for me. Thanks