count reads across exon junctions

[suninsky]

Suppose I have a bam file of paired-end reads and a gff file for exon locations. Is there an easy way to count the number of split reads (which can be identified by CIGAR string) and the number of paired-end reads that overlap any exon junctions? Assume we do not have any isoform information here...

Thanks a lot.

[ketan_bnf]

I am also finding a answer to this question, i am working on single end reads.

If somebody know about this please help.

[apfejes]

I think the question is not necessarily the right one. How do you know you've got your exon-junction-spanning reads spanning exons? I believe we usually align all of our transcriptome reads to reference genome + exon junction template, from which we can then reposition our exon-spanning reads.

If you haven't done that type of alignment, you may not even have your reads aligned correctly - but that would be incredibly dependent upon your aligner.

[ketan_bnf]

I think the question is not necessarily the right one. How do you know you've got your exon-junction-spanning reads spanning exons? I believe we usually align all of our transcriptome reads to reference genome + exon junction template, from which we can then reposition our exon-spanning reads.

If you haven't done that type of alignment, you may not even have your reads aligned correctly - but that would be incredibly dependent upon your aligner.

Sorry suninsky for posting in your thread,

I am having 454 exome reads aligned to ref seq using bwa to sam file.
Is it possible to find reads spanning exon region?

[apfejes]

Exome reads are usually aligned against the genome because they are captured from DNA, so, for the most part, you will not have exon-exon junctions in your data set. I'm not quite clear what it is you're asking.

You're probably best off trying to align as much as you cat to the genome and do your analysis on that set. If you're looking for reads spanning exon fusions, they'll probably be among the unaligned set, or if you want fusion proteins, you could use the paired end data to look for them.