Question: How to extract reads that map within a defined interval from each other
1
gravatar for dgimode
2.9 years ago by
dgimode60
dgimode60 wrote:

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

ADD COMMENTlink modified 2.9 years ago • written 2.9 years ago by dgimode60

Did you check bedtools closestBed

ADD REPLYlink written 2.9 years ago by venu6.2k
2
gravatar for Alex Reynolds
2.9 years ago by
Alex Reynolds28k
Seattle, WA USA
Alex Reynolds28k wrote:

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.

ADD COMMENTlink modified 2.9 years ago • written 2.9 years ago by Alex Reynolds28k
1
gravatar for dgimode
2.9 years ago by
dgimode60
dgimode60 wrote:

The BEDOPS way worked just fine for me. Thanks

ADD COMMENTlink written 2.9 years ago by dgimode60
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1664 users visited in the last hour