code: def density(args):
      
      p = OptionParser(density.__doc__)
      opts, args = p.parse_args(args)
  
      if len(args) != 2:
          sys.exit(not p.print_help())
  
      bedfile, fastafile = args
      bed = Bed(bedfile)
      sizes = Sizes(fastafile).mapping
      header = "seqid features size density_per_Mb".split()
      print("\t".join(header))
      for seqid, bb in bed.sub_beds():
          nfeats = len(bb)
          size = sizes[seqid]
          ds = nfeats * 1e6 / size
          print("\t".join(str(x) for x in \
                      (seqid, nfeats, size, "{0:.1f}".format(ds))))
docs: %prog density bedfile ref.fasta
  
      Calculates density of features per seqid.