How to get the output for each mapped read from TopHat?

[davidehs]

For the Tophat output file "junctions.bed", every line represents a single junction when mapping short sequencing reads to the genome. Each junction consists of two connected BED blocks, where each block is as long as the maximal overhang of any read spanning the junction (from TopHat manual).

For the following line,

chr20 257975 259040 JUNC00000002 6 - 257975 259040 255,0,0 2 49,70 0,995

it means 6 reads are mapped to this junction site. The limits of the two blockes, 257975 to 258024 (Block1) and 258970 to 259040 (Block2), were defined by any of the reads that spanning farthest at this junction site. Here is the problem. I want to know the start and end position for each of the 6 reads. So, the one line above may need to be expanded to 6 lines, with each line specifying the mapping position for each read.

The present Tophat can not generate the results I need. It seems that tophat produces some temporary files when running. Is there any way I can generate the results I need from these temporary file? Does anybody has similar experience, or knows any other program that can give the results I want?

Any suggestions will be greatly appreciated!

[ShaunMahony]

Parse the BAM output file for spliced reads that overlap your blocks.

[davidehs]

Hi ShaunMahony,

Thanks very much for your reply.

I am not familiar with BAM. Could you please give a bit more detailed explanation?

I appreciate it.


David

[swbarnes2]

I don't know TopHat, but BEDTools can take a bam file, and intersect it with a .bed file, and output the .bam entries that overlap the .bed file intervals.

Something like:

BEDTools/bin/intersectBed -abam TopHat.bam -b junctions.bed > junctions.bam

[ShaunMahony]

swbarnes2 -- yes that is a good idea, but the bed files produced by tophat describe splice junctions. So doing a simple overlap with a BAM could return all reads that are between the junction coordinates (including over introns and skipped exons), not just the reads that overlap the exon-exon boundary. I haven't used BEDTools, so I don't know how it will handle junction bed files.

davidehs -- I don't know of any other software tool to do exactly what you want, but it would be easy to write a simple script to parse the BAM file for junction reads that overlap block coordinates of interest (converting the BAM to SAM using samtools would allow you to text parse in perl or python, etc).

On the other hand, if you are just interested in a handful of junctions, why don't you just view the BAM file in a genome browser (UCSC?) and that should show you where all the reads start.

[davidehs]

It seems that I have to spend some time to get familiar with SAM, BAM, and BEDTools. They look pretty tough. Hope I can get through them soon.

Thanks for you guys help.

David