I don't know about tools, but i've used the following python code to calculate the ratio of reads that overlap the 5' or 3' ends of introns or that are spliced. We sum these across all introns in a gene set (we actaully use this for iCLIP analysis to see if RNA binding proteins bind pre-mRNA or spliced RNA).
import pysam
from collections import Counter
from CGAT import GTF, IOTools
def calculateSplicingIndex(bamfile, gtffile, outfile):
bamfile = pysam.AlignmentFile(bamfile)
counts = Counter()
for transcript in GTF.transcript_iterator(
GTF.iterator(IOTools.openFile(gtffile))):
introns = GTF.toIntronIntervals(transcript)
for intron in introns:
reads = bamfile.fetch(
reference=transcript[0].contig,
start=intron[0], end=intron[1])
for read in reads:
if 'N' in read.cigarstring:
blocks = read.get_blocks()
starts, ends = zip(*blocks)
if intron[0] in ends and intron[1] in starts:
counts["Exon_Exon"] += 1
else:
counts["spliced_uncounted"] += 1
elif (read.reference_start <= intron[0] - 3
and read.reference_end >= intron[0] + 3):
if transcript[0].strand == "+":
counts["Exon_Intron"] += 1
else:
counts["Intron_Exon"] += 1
elif (read.reference_start <= intron[1] - 3
and read.reference_end >= intron[1] + 3):
if transcript[0].strand == "+":
counts["Intron_Exon"] += 1
else:
counts["Exon_Intron"] += 1
else:
counts["unspliced_uncounted"] += 1
header = ["Exon_Exon",
"Exon_Intron",
"Intron_Exon",
"spliced_uncounted",
"unspliced_uncounted"]
with IOTools.openFile(outfile, "w") as outf:
outf.write("\t".join(header)+"\n")
outf.write("\t".join(map(str, [counts[col] for col in header]))
+ "\n")
Unfortunately this uses a bunch of libraries you may or may not have, including CGAT (for the GTF parser and IOTools package) and pysam.
Once you've got these statistics you can calculate the "splicing index" as
the log2 ratio of 2 times the number of spliced reads divided by the number of reads overlapping the 3' and 5' ends of introns.