Finding reads shared between 2 bam files

[plumb_r]

Hi
I have mapped bisulfite converted HiSeq reads separately to in silico converted forward & reverse strands so have files

sample.forward.bam
sample.reverse.bam

I now would like a count of the reads which map to both and only both of the strands i.e. those which are in both forward & reverse bam files but not in only one of them.
Any ideas?
I was thinking of using either:

samtools merge out.bam sample.forward.bam sample.reverse.bam

but can't work out how this would help
Could another option could be to use bedtools to convert both bam to bed & then intersect these? Something like:

bamToBed -i sample.forward.bam > for.reads.bed
bamToBed -i sample.reverse.bam > rev.reads.bed
intersectBed -a .for.reads.bed -b rev.reads.bed -c

Would this work?
Thanks

[dpryan]

This will somewhat depend on how you did the alignments. You also run into the problem of whether to consider a read that just barely mapped to say, the original top strand (assume a large edit distance), but mapped perfectly with no edits to the original bottom strand as actually mapping to both strands or not.

Converting the alignments to BED files and then intersecting them won't do what you want, since you'll be intersecting covered areas rather than whether an individual read maps (e.g. if a read aligns to chr1 of the top strand and chr2 of the bottom strand then the intersection (which is coordinate based) won't work). You could merge the alignments and name sort them and then just write a program to process that. That would be quite straight forward. You could also use an almost identical program on the unmerged files to save a step (and some time).

Might I ask what you're trying to accomplish with this? If you're trying to just figure out which of your reads to exclude from downstream analysis because of ambiguous mapping then you could just rerun the alignments with Bismark or Bison (of which I'm the author).