Given an bam file, find all potential introns, optionally within a given region.

countJuctions.py

#!/usr/bin/env python
import pysam
import argparse
import collections

parser = argparse.ArgumentParser()
parser.add_argument('bamfile',
                    help="bamfile for analysis")
parser.add_argument('-r','--region',
                    help="Limit to region chrN:XXX-YYY")

opts=parser.parse_args()

samfile = pysam.Samfile(opts.bamfile)
reads = samfile.fetch()
print opts.region
if(opts.region):
    reads = samfile.fetch(region=opts.region)

introns=collections.defaultdict(int)
for read in reads:
    currentloc = read.pos
    for cigarop in read.cigar:
        if(cigarop[0]==3):
            introns[(samfile.getrname(read.tid),currentloc,currentloc+cigarop[1])]+=1
        currentloc=currentloc+cigarop[1]

for i in introns:
    print i,introns[i]