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
Sure, you can do this with a few genomic set operations with [BEDOPS][1] 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.

  [1]: https://github.com/bedops/bedops