This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extract reads overlapping a specific region in bam file

Hello,

ref --------------start-----------------------------------------stop-------------------------------
r1                   -------------------------------------------------------
r2             ---------------------------------------------------------------
r3                  -----------------    --------------------------------------
r4              -----------------------------------------------------------------------------

a. I would like to extract reads from bam file that overlap entirely the start and stop region from both directions (from start to stop and from stop to start). From the example above, I need to keep only r1,r2 and r4, but not r3.

How to do this?

I tried this, but I still have some small fragmented reads...

samtools view -b -h -q 10 input.bam chrX:230-330 | awk 'BEGIN{OFS="\t"}{if($1 ~ /^"@"/) {print} else {if($4 >= 230) {print} else {}}}' | samtools view -Sbo output.bam

I also tried:

samtools view -h -q 10 input.bam chrX:230-330 | awk 'BEGIN{OFS="\t"}{if($1 ~ /^"@"/) {print} else {if($4 >= 230 && length($10) >= 100) {print} else {}}}' | samtools view -Sbo output.bam -
awk samtools
• 4,963 views • link
updated by hopf • hopf • written by User000 User000

Your diagram and text do not match - r1 is the only read that spans both the start and stop. Is that what you want?

• 0 views • link
written by harold.smith.tarheel harold.smith.tarheel

I updated my question, thanks

• 0 views • link
written by User000 User000

If the regions are small you could combine the standard region querying with the filter expression, ie

samtools view -e 'pos <= 230 && pos+rlen >= 330' in.bam chrX:230-330

I haven't tested it. Maybe it needs to be pos+rlen-1, but you get the jist. I hope it works! If not then maybe it'll need a samtools view | samtools view style command to do it in two passes.

• 0 views • link
written by jkbonfield jkbonfield

1 answer

Create a BED file that contains the start and stop positions, then use BEDtools:

bedtools intersect -wa BAM -b BED -F 1.0
• 0 views • link
written by harold.smith.tarheel harold.smith.tarheel

The user wants reads that 100% / completely fall within a desired interval. Does this BEDTools command do that? I guess that, yes, it can do it, and that -F 1.0 is the key.

-f Minimum overlap required as a fraction of A. Default is 1E-9 (i.e. 1bp).

-F Minimum overlap required as a fraction of B. Default is 1E-9 (i.e., 1bp).

User000, can you please try this?

• 0 views • link
written by Kevin Blighe Kevin Blighe

I've tried this, and I get the error:

***** ERROR: Unrecognized parameter: input.bam *****

Also if i try bedtools intersect -wa -a BAM -b BED -F 1.0, I get this error:

ERROR: Received illegal bin number 4294967295 from getBin call.
ERROR: Unable to add record to tree.

Also, I have only one interval of interest, not a multiple, so my bed file contatins only that. It that ok?

• 0 views • link
written by User000 User000

What are the names of your BED and BAM files?

• 0 views • link
written by Kevin Blighe Kevin Blighe

Try:

bedtools intersect -wa -a BAM -b BED -F 1.0
• 0 views • link
written by hopf • hopf

Log in to answer this question.

More posts like this