How to get reads for a particular splice junction – TopHat and BEDTools

[ggould]

Does anyone have a quick way of getting all the reads that correspond to a particular splice junction? I’ve mapped my RNA-seq data with TopHat, made a new .bam file with only reads that span junctions (N in the cigar string), and run intersectBed with this new junction-only bam file and the junctions.bed file created by TopHat. The issue I’m encountering is if an intron has an alternative splice site, say this particular intron has only 5’ splice site and two alternative 3’ splice sites, a read will “intersect” this junction if it crosses the 5’ss but that read might not belong to that junction because it uses a different 3’ splice site. Is there a quick way of pulling out only the reads that span a particular junction using BEDTools (or another program) or am I better off doing this myself in python? Thanks!

[GenoMax]

How about the "samtools view region" option?

Usage:   samtools view [options] <in.bam>|<in.sam> [region1 [...]]

Options: -b       output BAM
         -h       print header for the SAM output
         -H       print header only (no alignments)
         -S       input is SAM
         -u       uncompressed BAM output (force -b)
         -1       fast compression (force -b)
         -x       output FLAG in HEX (samtools-C specific)
         -X       output FLAG in string (samtools-C specific)
         -c       print only the count of matching records
         -B       collapse the backward CIGAR operation
         -@ INT   number of BAM compression threads [0]
         -L FILE  output alignments overlapping the input BED FILE [null]
         -t FILE  list of reference names and lengths (force -S) [null]
         -T FILE  reference sequence file (force -S) [null]
         -o FILE  output file name [stdout]
         -R FILE  list of read groups to be outputted [null]
         -f INT   required flag, 0 for unset [0]
         -F INT   filtering flag, 0 for unset [0]
         -q INT   minimum mapping quality [0]
         -l STR   only output reads in library STR [null]
         -r STR   only output reads in read group STR [null]
         -s FLOAT fraction of templates to subsample; integer part as seed [-1]
         -?       longer help

[sdriscoll]

what if you focused on the two exons on either side of the intron rather than the intron itself? certainly if a spliced alignment hits BOTH of those exons you've got a winner. in fact if you even know the coordinates of the intron you could just fake a bed file with two rows with 10 or 20bp on either side of the intron and then intersect that with your BAM using 'bedtools intersect' then parse for cases were a single read hit both of them.