I am wondering if it is possible via samtools/bamtools to extract all the paired-end mapped reads (against the human hg19 genome) that fall in exon regions. I have a BED file that defines all these regions. I already know that it's possible to calculate this ($ coverageBed -abam reads.bam -b exons.bed > exons.bed.coverage), but I have no clue how to extract the reads from the .bam file). I tried using the samtools/bamtools manual but couldn't figure it out.
Anyone who can help me?