[BioC] Extract junction reads from bam

[Helen Zhou]

Dear Sir/Madam,

I have a bam file with reads from a paired-end RNA-seq experiments, and I would like to extract all reads that span a particular intron/exon junction. For example where part of a read maps to exon ABC starting at chr 1 position 66909348. I do not know where the first part of the read maps; there are multiple possible exons being spliced to exon ABC.

I can extract all reads mapping to and spanning for example 50nt around the junction, saying:
library(GenomicRanges)
loc <-RangesList('1'=IRanges(start=6690298,end=6690398)) 
parameter <- ScanBamParam(which=loc, what=extract, simpleCigar=FALSE, reverseComplement=FALSE)
aln <- readGAlignmentsFromBam("test.bam", param=parameter)

However, this includes all the reads that span across this region because the ends of the read map outside the area indicated in 'loc'. How can I get only the reads where (part of) the read starts mapping at the exact intron/exon location at position 66909348?

Thanks you
H Zhou
	[[alternative HTML version deleted]]