code: def insertion(args):
      
      p = OptionParser(insertion.__doc__)
      p.add_option("--mindepth", default=6, type="int",
                   help="Minimum depth to call an insertion")
      p.set_outfile()
      opts, args = p.parse_args(args)
  
      if len(args) != 1:
          sys.exit(not p.print_help())
  
      bedfile, = args
      mindepth = opts.mindepth
      bed = Bed(bedfile)
      fw = must_open(opts.outfile, "w")
      for seqid, feats in bed.sub_beds():
          left_ends = Counter([x.start for x in feats])
          right_ends = Counter([x.end for x in feats])
          selected = []
          for le, count in left_ends.items():
              if count >= mindepth:
                  selected.append((seqid, le, "LE-{0}".format(le), count))
          for re, count in right_ends.items():
              if count >= mindepth:
                  selected.append((seqid, re, "RE-{0}".format(re), count))
          selected.sort()
          for seqid, pos, label, count in selected:
              label = "{0}-r{1}".format(label, count)
              print("\t".join((seqid, str(pos - 1), str(pos), label)), file=fw)
docs: %prog insertion mic.mac.bed
  
      Find IES based on mapping MIC reads to MAC genome. Output a bedfile with
      'lesions' (stack of broken reads) in the MAC genome.