code: def nucmer(args):
      
      p = OptionParser(nucmer.__doc__)
      opts, args = p.parse_args(args)
  
      if len(args) != 5:
          sys.exit(not p.print_help())
  
      mapbed, mtrfasta, asmfasta, chr, idx = args
      idx = int(idx)
      m1 = 1000000
      bedfile = "sample.bed"
      bed = Bed()
      bed.add("\t".join(str(x) for x in (chr, (idx - 1) * m1, idx * m1)))
      bed.print_to_file(bedfile)
  
      cmd = "intersectBed -a {0} -b {1} -nonamecheck -sorted | cut -f4".\
                  format(mapbed, bedfile)
      idsfile = "query.ids"
      sh(cmd, outfile=idsfile)
  
      sfasta = fastaFromBed(bedfile, mtrfasta)
      qfasta = "query.fasta"
      cmd = "faSomeRecords {0} {1} {2}".format(asmfasta, idsfile, qfasta)
      sh(cmd)
  
      cmd = "nucmer {0} {1}".format(sfasta, qfasta)
      sh(cmd)
  
      mummerplot_main(["out.delta", "--refcov=0"])
      sh("mv out.pdf {0}.{1}.pdf".format(chr, idx))
docs: %prog nucmer mappings.bed MTR.fasta assembly.fasta chr1 3
  
      Select specific chromosome region based on MTR mapping. The above command
      will extract chr1:2,000,001-3,000,000.