Count Junctions In A Sam/Bam File

Question

[Dave Bridges]

When considering transcript variants, I like using DEXseq to count reads mapped to exons, but while that is useful, I thought it would be more useful to be able to count junctions.

For example a given gene structure could be

1-2-3-4-5

But there might also be these structures

1-2-4-5

1-3-4-5

Although DEXseq seems like it could provide insight on this, in the above example, if both transcripts are about the same expression DEXseq would give similar counts for exons 2 and 3 and it would be hard to prove that these transcripts exist.

What I think would be better is if there would be a way to count junctions in a BAM file, such that you could say something like:

80% of transcripts have junctions from 1-2

20% of transcripts have junctions from 1-3

For a given junction. I think cufflinks tries to do something like this but it is not clear how they get their whole transcript numbers from these small data sets. Currently, for genes of interest I am loading the bam files into IGV and counting exon junctions by hand, but is there a script or tool that does this, maybe only counting junction transcripts with reasonable overlaps to both sides of the junction?

Answer 1

[Malachi Griffith]

If you ran TopHat to produce your BAM file it will have produced a 'junctions.bed' file. The score column of this junctions file contains the read counts for each junction. Since it is a BED file it can be manipulated with existing tools such as BEDTools. It can also be loaded in IGV alongside your BAM file.

By extracting junction coordinates from this file you can easily do things like:

• Determining which junctions involve canonical GT-AG splice sites or other splice sites
• Determine intron size
• Identify junctions that correspond to known transcripts, novel junctions that are anchored to a known acceptor or donor site, novel junctions that involve known acceptor and donor sites, but in a novel combination (i.e. exon skipping)
• For exon skipping events determine the number of exons skipped
• etc.

Answer 2

[Dan D]

I think FastMISO is what you want. It was a bear to install for me, but I've been happy with it thus far.

Answer 3

[Ryan Thompson]

You might want to try the spliceGraph function in the GenomicFeatures package of BioConductor. Here is a PDF files that gives a detailed example of how to use it with DEXSeq:

http://www.bioconductor.org/packages//2.10/bioc/vignettes/GenomicFeatures/inst/doc/spliceGraph.pdf

I haven't tried it myself yet, but I plan to.

Answer 4

[madk00k]

New version of Qualimap tool counts junctions in RNA-seq Quality Control mode.

QualiMap can downloaded from here.

Here's an example of report.