code: def juncs(args):
      
      from tempfile import mkstemp
      from pybedtools import BedTool
  
      p = OptionParser(juncs.__doc__)
      p.set_outfile()
  
      opts, args = p.parse_args(args)
  
      if len(args) < 1:
          sys.exit(not p.print_help())
  
      fh, trimbed = mkstemp(suffix = ".bed")
      fw = must_open(trimbed, "w")
      for i, juncbed in enumerate(args):
          bed = Bed(juncbed, juncs=True)
          for b in bed:
              ovh = [int(x) for x in b.extra[-2].split(",")]
              b.start += ovh[0]
              b.end -= ovh[1]
              b.accn = "{0}-{1}".format(b.accn, i)
              b.extra = None
              print(b, file=fw)
      fw.close()
  
      if len(args) > 1:
          sh("sort -k1,1 -k2,2n {0} -o {0}".format(trimbed))
  
          tbed = BedTool(trimbed)
          grouptbed = tbed.groupby(g=[1,2,3,6], c=5, ops=['sum'])
  
          cmd = """awk -F $'\t' 'BEGIN { OFS = FS } { ID = sprintf("mJUNC%07d", NR); print $1,$2,$3,ID,$5,$4; }'"""
          infile = grouptbed.fn
          sh(cmd, infile=infile, outfile=opts.outfile)
      else:
          sort([trimbed, "-o", opts.outfile])
  
      os.unlink(trimbed)
docs: %prog junctions junctions1.bed [junctions2.bed ...]
  
      Given a TopHat junctions.bed file, trim the read overhang to get intron span
  
      If more than one junction bed file is provided, uniq the junctions and
      calculate cumulative (sum) junction support