Finding Reads Joining Exons From Tophat Output

Question

[hbw]

I did an alignment with tophat, supplying the GFF file for gene annotations, and the chromosomes fasta (in bowtie2 index format). I converted the output accepted_hits.bam to SAM. In the SAM header, it is listing all the chromosomes, and the reads are mapped to the chromosmes, like this:

266GRCF1GCCAAT:1:1101:8591:43171        0 NC_000076.6   38894021        50      49M     *       0       0       GGCGGGGGGAGGGTATCATGGACTTTTGGGATAGCATTTGAAATATAAA      @@@DDDDDB6BBB+8>CDDCCCCCCCCCBB?ACCCCECCCCCCACCCDE AS:i:0 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:49 YT:Z:UU NH:i:1
266GRCF1GCCAAT:1:2211:7290:90876        0 NC_000076.6   38894073        50      49M     *       0       0       AGAAAATATGTTAAAAATGTTTTGGAAAAAAGAAAAAAAGAAAATTGCC      CCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJJJJJIJJJJJJJJJJJ AS:i:0 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:49 YT:Z:UU NH:i:1
266GRCF1GCCAAT:1:2202:7215:55833        0 NC_000076.6   38894074        50      49M     *       0       0       GAAAATATGTTAAAAATGTTTTGGAAAAAAGAAAAAAAGAAAATTGCCA      @CCFFFFFHHHHHJJJJIGIJJJIJJIIIIGIIIIGGEGHJIJIEGHII AS:i:0 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:49 YT:Z:UU NH:i:1

Where are the read alignments that map to two exons separated by an intron? How do I extract this information from the Tophat output?

Answer 1

[Sean Davis]

You are looking for reads that have the "N" operation in the samtools CIGAR column. See the SAM specification for details.

Answer 2

[Charles Warden]

Have you checked the junctions.bed output file?

I think this should contain the sort of statistics that you are looking for (counts for reads supporting a given exon junction)