Samtools or Bedtools: How to filter a bam file with a bed file using strand information

Question

[retrogenomics]

Hi

I would like to filter a bam file, keeping only reads overlapping with genomic intervals from a bed file. I used samtools for this:

samtools view -b -h -L bedfile.bed bamfile.bam

However the -L option does not seem to take into account the strand information.

Do you know if there is another option or way to do it that would keep strand information?

[Devon Ryan]

Do all of the regions in the BED file have the same orientation? If not, you'll be better off with bedtools or bedops.

[retrogenomics]

Some intervals in the bed files are +, others are -. With bed tools, I don't see an easy way to manipulate the bam file. The output I'm looking for is an alignment that I can visualize in IGV.

[Devon Ryan]

bedtools intersect will output a BAM file if you give it one as input. Just sort and index it and then load into IGV.

Answer 1

[swbarnes2]

The kind of brute force way to do it is to split apart the + and- strand intervals, and use samtools view to filter the .bam by direction before piping into BEDTools, then merge the two resulting files.

Answer 2

[Alex Reynolds]

Another option:

$ samtools view -F 0x10 reads.bam | bam2bed > reads.forward.bed
$ samtools view -f 0x10 reads.bam | bam2bed > reads.reverse.bed
$ bedops --everything reads.*.bed > reads.bed

Then filter or map reads.bed with bedops or bedmap, as usual, to generate subsets or new sets in BED format. IGV can support BED files.