<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def simulate(args):
 """ %prog simulate run_dir 1 300 Simulate BAMs with varying inserts with dwgsim. The above command will simulate between 1 to 300 CAGs in the HD region, in a directory called `run_dir`. """

    p = OptionParser(simulate.__doc__)
    p.add_option("--method", choices=("wgsim", "eagle"), default="eagle",
                 help="Read simulator")
    p.add_option("--ref", default="hg38", choices=("hg38", "hg19"),
                 help="Reference genome version")
    p.add_option("--tred", default="HD", help="TRED locus")
    add_simulate_options(p)
    opts, args = p.parse_args(args)

    if len(args) != 3:
        sys.exit(not p.print_help())

    rundir, startunits, endunits = args
    ref = opts.ref
    ref_fasta = "/mnt/ref/{}.upper.fa".format(ref)
    startunits, endunits = int(startunits), int(endunits)
    basecwd = os.getcwd()
    mkdir(rundir)
    os.chdir(rundir)
    cwd = os.getcwd()

    # TRED region (e.g. Huntington)
    pad_left, pad_right = 1000, 10000
    repo = TREDsRepo(ref=ref)
    tred = repo[opts.tred]
    chr, start, end = tred.chr, tred.repeat_start, tred.repeat_end

    logging.debug("Simulating {}".format(tred))
    fasta = Fasta(ref_fasta)
    seq_left = fasta[chr][start - pad_left:start - 1]
    seq_right = fasta[chr][end: end + pad_right]
    motif = tred.repeat

    simulate_method = wgsim if opts.method == "wgsim" else eagle
    # Write fake sequence
    for units in range(startunits, endunits + 1):
        pf = str(units)
        mkdir(pf)
        os.chdir(pf)
        seq = str(seq_left) + motif * units + str(seq_right)
        fastafile = pf + ".fasta"
        make_fasta(seq, fastafile, id=chr.upper())

        # Simulate reads on it
        simulate_method([fastafile, "--depth={}".format(opts.depth),
                          "--readlen={}".format(opts.readlen),
                          "--distance={}".format(opts.distance),
                          "--outfile={}".format(pf)])

        read1 = pf + ".bwa.read1.fastq"
        read2 = pf + ".bwa.read2.fastq"
        samfile, _ = align([ref_fasta, read1, read2])
        indexed_samfile = index([samfile])

        sh("mv {} ../{}.bam".format(indexed_samfile, pf))
        sh("mv {}.bai ../{}.bam.bai".format(indexed_samfile, pf))

        os.chdir(cwd)
        shutil.rmtree(pf)

    os.chdir(basecwd)